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ABSTRACT 

We simulate the multi-phase interstellar medium (ISM) randomly heated and stirred by su- 
pemovae (SNe), with gravity, differential rotation and other parameters of the solar neigh- 
bourhood. Here we describe in detail both numerical and physical aspects of the model, in- 
eluding injection of thermal and kinetic energy by SN explosions, radiative cooling, photo- 
electric heating and various transport processes. With a three-dimensional domain extending 
4^ 1x1 kpc^ horizontally and 2 kpc vertically (symmetric about the galactic mid-plane), the 

Qh model routinely spans gas number densities 10~^-10^ cm~^, temperatures 10-10^ K, local 

Q velocities up to 10^ km (with Mach number up to 25). The working numerical resolution 

^ of 4 pc has been selected via simulations of a single expanding SN remnant, where we closely 

^ reproduce, at this resolution, analytical solutions for the adiabatic and snowplough regimes. 

The feedback of the the halo on the disc cannot be captured in our model where the domain 
only extends to the height of 1 kpc above the mid-plane. We argue that to reliably model the 
disc-halo connections would require extending the domain horizontally as well as vertically 
due to the increasing horizontal scale of the gas flows with height. 

The thermal structure of the modelled ISM is classified by inspection of the joint probabil- 
ity density of the gas number density and temperature. We confirm that most of the complexity 
V/^ can be captured in terms of just three phases, separated by temperature borderlines at about 

^T) 10^ K and 5 x 10^ K. The probability distribution of gas density within each phase is approx- 

imately lognormal. We clarify the connection between the fractional volume of a phase and 
its various proxies, and derive an exact relation between the fractional volume and the filling 
^vq factors defined in terms of the volume and probabilistic averages. These results are discussed 

^-H in both observational and computational contexts. The correlation scale of the random flows 

r* is calculated from the velocity autocorrelation function; it is of order 100 pc and tends to grow 

. with distance from the mid-plane. We use two distinct parameterizations of radiative cooling 

to show that the multi-phase structure of the gas is robust, as it does not depend significantly 
^ on this choice. 
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1 INTRODUCTION 

The multi-phase structure of the interstellar medium (ISM) affects 
almost all aspects of its dynamics, including its evolution, star for- 
mation, galactic winds and fountains, and the behaviour of mag- 
netic fields and cosmic rays. In a widely accepted picture (Cox & 
Smith 1974; McKee & Ostriker 1977), most of the volume is oc- 
cupied by the hot (T 10^ K), warm (T 10^ K) and cold 
(T c::^ 10^ K) phases. The concept of the multi-phase ISM in pres- 
sure equilibrium has endured with modest refinement (Cox 2005), 
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e.g., deviations from thermal pressure balance have been detected 
(Kalberla & Kerp 2009, and references therein). Dense molecular 
clouds, while binding most of the total mass of the interstellar gas 
and being of key importance for star formation, occupy a negligi- 
ble fraction of the total volume (e.g. Kulkarni & Heiles 1987, 1988; 
Spitzer 1990; McKee 1995). The main sources of energy maintain- 
ing this complex structure are supernova (SN) explosions and stel- 
lar winds (Mac Low & Klessen 2004, and references therein). The 
clustering of SNe in OB associations facilitates the escape of the 
hot gas into the halo thus reducing the volume filling factor of the 
hot gas in the disc, perhaps down to 10% at the mid-plane (Norman 
& Ikeuchi 1989). The energy injected by the SNe not only produces 
the hot gas but also drives ubiquitous compressible turbulence in all 
phases, as well as driving outflows from the disc, associated with 
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the galactic fountain or wind, as first suggested by Bregman (1980). 
Thus turbulence, the multi-phase structure, and the disc-halo con- 
nection are intrinsically related features of the ISM. 

A comprehensive description of the complex dynamics of 
the multi-phase ISM has been significantly advanced by numer- 
ical simulations in the last three decades, starting with Chiang 
& Prendergast (1985), followed by many others including Rosen 
et al. (1993); Rosen & Bregman (1995); Vazquez- Semadeni et al. 
(1995); Passot et al. (1995); Rosen et al. (1996); Korpi et al. (1999); 
Gazol-Patino & Passot (1999); Wada & Norman (1999); de Avillez 
(2000); Wada & Norman (2001); de Avillez & Berry (2001); de 
Avillez & Mac Low (2002); Wada et al. (2002); de Avillez & Bre- 
itschwerdt (2004); Balsara et al. (2004); de Avillez & Breitschwerdt 
(2005a);de Avillez & Breitschwerdt (2005b); Slyz et al. (2005); 
Mac Low et al. (2005); Joung & Mac Low (2006); de Avillez 
& Breitschwerdt (2007); Wada & Norman (2007); Gressel et al. 
(2008). Numerical simulations of this type are demanding even 
with the best computers and numerical methods available. The self- 
regulation cycle of the ISM includes physical processes spanning 
enormous ranges of gas temperature and density, and of spatial 
and temporal scales, as it involves star formation in the cores of 
molecular clouds, assisted by gravitational and thermal instabili- 
ties at larger scales, which evolve against the global background of 
transonic turbulence driven, in turn, by star formation (Mac Low 
& Klessen 2004). It is understandable that none of the existing nu- 
merical models covers the whole range of parameters, scales and 
physical processes known to be important. 

Two major approaches in earlier work focus either on the dy- 
namics of diffuse gas or on dense molecular clouds. Our model 
belongs to the former class, where we are mainly concerned with 
the ISM dynamics in the range of scales of order 10 pc-1 kpc. Nu- 
merical constraints prevent us (like many other authors) from fully 
including the gravitational and thermal instabilities which involve 
scales of less than 1 pc. In order to assess the sensitivity of our re- 
sults to the parameterization of radiative coohng, we consider mod- 
els with thermal instability, but reduce its efficiency using a suffi- 
ciently strong thermal conductivity to avoid the emergence of struc- 
tures that are unresolvable at our numerical resolution. The results 
are compared to models with no thermally unstable branch over 
the temperature range between the cold and warm phases. To our 
knowledge, no direct study addressing the difference between these 
two kinds of the cooling parameterizations has been made. We note, 
however, that Vazquez-Semadeni et al. (2000) compared their ther- 
mally unstable model to a different model by Scalo et al. (1998), 
who used a thermally stable cooling function. Similarly, de Avillez 
& Breitschwerdt (2004) and Joung & Mac Low (2006) compared 
results obtained with different cooling functions, but again com- 
paring different models: here we compare models with different 
cooling functions but which are otherwise the same. 

An unavoidable consequence of the modest numerical resolu- 
tion available, if we are to capture the dynamics on 1 kpc scales, 
is that star formation, manifesting itself only through the ongoing 
SN activity in our model, has to be heavily parameterized. We do, 
however, ensure that individual supernova remnants are modelled 
accurately, since this is essential to reliably reproduce the injection 
of thermal and kinetic energy into the ISM. In particular, our model 
reproduces with high accuracy the evolution of supernova remnants 
from the Sedov-Taylor stage until the remnant disintegrates and 
merges into the ISM (Appendix B). 

The dimensionless parameters characteristic of the ISM, such 
as the kinetic and magnetic Reynolds numbers (reflecting the rel- 
ative importance of gas viscosity and electrical resistivity) and the 



Prandtl number (quantifying thermal conductivity), are too large to 
be simulated with current computers. Similarly to most numerical 
simulations of this complexity, our numerical techniques involve a 
range of artificial transport coefficients for momentum and thermal 
energy (such as shock-capturing viscosity). We explore and report 
here the sensitivity of our results to the artificial elements in our 
basic equations. 

This paper is the first of a planned series, in which we aim to 
clarify which components and physical processes control the dif- 
ferent properties of the ISM. Our next step is to add magnetic fields 
to the model, to study both their origin and role in shaping the ISM. 
But in order to identify where the magnetic field is important and 
where it is not, we first must understand what the properties of a 
purely hydrodynamic ISM would be. 

The structure of the paper is as follows. In Section 2 we 
present our basic equations, numerical methods, initial and bound- 
ary conditions, as well as the physical ingredients of the model, 
such as our modelling of SN activity and heating and cooling of 
the ISM. Our results are presented in Sections 3-8, including an 
overview of the multi-phase structure of the ISM, the correlation 
length of random flows, and their sensitivity to the cooling func- 
tion and numerical resolution. Our results are discussed in a broader 
context in Section 9, where our conclusions are also summarised. 
Detailed discussion of important technical and numerical aspects of 
the model, and the effects of the unavoidable unphysical assump- 
tions adopted, can be found in Appendices: the accuracy of our 
modelling of individual supernova remnants in Appendix B, our 
control of numerical dissipation in Appendix C, and sensitivity to 
thermal instability in Appendix D. 



2 BASIC EQUATIONS AND THEIR NUMERICAL 
IMPLEMENTATION 

2.1 Basic equations 

We solve numerically a system of hydrodynamic equations using 
the Pencil Code (http://code.google.eom/p/pencil-code) which is 
designed for fully nonlinear, compressible magnetohydrodynamic 
(MHD) simulations. We consider only the hydrodynamic regime 
for the purposes of this paper; MHD simulations, which are in 
progress, will be reported elsewhere. Nor do we include cosmic 
rays, which we subsequently plan to add to the MHD simulations. 

The basic equations include the mass conservation equation, 
the Navier-Stokes equation (written here in the rotating frame), and 
the heat equation written in terms of the specific entropy:^ 

^ = -pV-tx + psN, (1) 

Du _i_ 2^ / / IN 

- -p VcrsN - Cs V (s/cp + Inp) 

- V$ - Su^y -2Q.XU 

+ i/(V^tx+ |VV-tt + 2W. Vlnp) 

+ C.(VV.^), (2) 

pT^ = asN + pT - p'A + V • (cppxVT) + 2pv |W|' 

+ CxP(V.^)% (3) 

^ For the reader's convenience, Appendix A contains a list of variables 
used in the text with their definitions. 
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where p, T and s are the gas density, temperature and specific en- 
tropy, respectively, u is the deviation of the gas velocity from the 
background rotation profile (here called the velocity perturbation), 
Cs is the adiabatic speed of sound, Cp is the heat capacity at constant 
pressure, S is the velocity shear rate associated with the Galactic 
differential rotation at the angular velocity ft assumed to be aligned 
with the 2:-axis (see below). The Navier-Stokes equation includes 
viscosity u and the rate of strain tensor W whose components are 
given by 

dui du 



dxj dxi 3 ' 



(4) 



as well as the shock-capturing viscosity ^i^. The system is driven 
by SN energy injection, at the rates &sn (per unit volume) in the 
form of kinetic energy in Eq. (2) and thermal energy in Eq. (3). 
Energy injection is applied in a single time step and is confined to 
the interiors of newly introduced SN remnants, and the total energy 
injected per supernova is denoted Eq^. The mass of the SN ejecta 
is included in Eq. (1) via the source yOsN- The forms of these terms 
are specified and further details are given in Section 2.2. The heat 
equation also contains a thermal energy source due to photoelectric 
heating pF, energy loss due to optically thin radiative cooling A, 
heat conduction with the thermal diffusivity x (with K = Cppx the 
thermal conductivity), viscous heating (with |W| the determinant 
of W), and the shock-capturing thermal diffusivity Cx- 
The advective derivative, 

§i = l + iU + u).V, (5) 

includes transport by an imposed shear flow U = (0, Sx, 0) in 
the local Cartesian coordinates (taken to be linear across the local 
simulation box), with the velocity u representing a deviation from 
the overall rotational velocity U. As will be discussed later, due 
to anisotropies (e.g. density stratification, anisotropic turbulence), 
large-scale flows will be generated in the system; one example is 
the systematic vertical outflow discussed at length in this paper. 
Therefore, the perturbation velocity u consists of two parts, a mean 
flow and random velocities. Here we consider a mean flow obtained 
by Gaussian smoothing (Germano 1992): 



{u)^{x) = / u{x)Gi{x — x)d^x, 
Jv 



(6) 



where we use a smoothing scale i ^ 50 pc, necessarily some- 
what shorter than the flow correlation length Zq obtained in Sec- 
tion 6 (for details, see Gent et al. 2013). The random flow is then 
ifco = — a- The differential rotation of the galaxy is mod- 
elled with a background shear flow along the local azimuthal {y) 
direction, JJy — Sx. The shear rate is S = r9Q/9r in terms of 
galactocentric distance r, which translates into the x -coordinate 
of the local Cartesian frame. In this paper we consider models 
with rotation and shear similar to those in the solar neighbourhood, 
— —S = 25kms~^ kpc~^. We do not expect the gas veloci- 
ties and thermal structure discussed here to depend strongly on the 
rotation and shear parameters, although other aspects of the solu- 
tion will be more sensitive to these. Future papers will consider the 
rotation and shear in more detail; and will also include magnetic 
fields, whose generation may depend strongly on these parameters. 
We consider an ideal gas, with thermal pressure given by 

P = pT, 

pnip 

where /cb is the Boltzmann constant, nip is the proton mass, and 



p — 0.62 is the mean molecular weight of a fully ionised gas of 
the Solar chemical composition. 

In Eq. (2), $ is the gravitational potential produced by stars 
and dark matter. For the Solar vicinity of the Milky Way, Kuijken 
& Gilmore (1989) suggest the following form of the vertical gravi- 
tational acceleration (see also Ferriere 2001): 



9$ 

9z = -w- = - 



ai 



dz 



(7) 



withai = 4.4 x 10"^^kms-^a2 = 1.7 x 10"^^ kms-^ = 
200 pc and Z2 = 1 kpc. We neglect self-gravity of the interstellar 
gas because it is subdominant at the scales of interest. 



2.2 Modelling supernova activity 

We include both Type II and Type I SNe in our simulations, dis- 
tinguished only by their frequency and vertical distribution. The 
SNe frequencies are those in the Solar neighbourhood (e.g. Tam- 
mann, Loffler & Schroder 1994). Type II SNe are introduced at a 
rate, per unit surface area, of vu — 25kpc~^ Myr~^ (0.02 yr~^ 
in the whole Galaxy), with fluctuations of the order of 10~^ yr~^ 
at a time scale of order 10 Myr. Such fluctuations in the Type II 
SN rate are natural to introduce; there is some evidence that they 
can enhance dynamo action in MHD models (Hanasz et al. 2004; 
Balsara et al. 2004). The surface density rate of Type I SNe is 
h>i — 4kpc~^ Myr~^ (interval of 290 years between Type I SN 
explosions in the Galaxy). We do not explicitly include any spatial 
clustering of the SNe. 

Unlike most other ISM models of this type, the SN energy in 
the injection site is split between thermal and kinetic parts, in or- 
der to reduce artificial temperature and energy losses at early stages 
of the SN remnant evolution. Thermal energy density is distributed 
within the injection site as exp[— (r/rsN)^], with r the local spher- 
ical radius and rsN (of order 10 pc - see below) the nominal loca- 
tion of the remnant shell (i.e. the radius of the SN bubble) at the 
time of injection. Kinetic energy is injected by adding a spherically 
symmetric velocity field Ur oc exp[— (r/rsN)^]; subsequently, this 
rapidly redistributes matter into a shell. To avoid a discontinuity in 
u at the centre of the injection site, the centre is simply placed mid- 
way between grid points. We also inject 4M0 as stellar ejecta, with 
density profile exp[— (r/rsN)^]. Given the turbulent environment, 
there are significant random motions and density inhomogeneities 
within the injection regions. Thus, the initial kinetic energy is not 
the same in each region, and, injecting part of the SN energy in the 
kinetic form results in the total kinetic energy varying between SN 
remnants. We therefore record the energy added for every remnant 
so we can fully account with the rate of energy injection. For ex- 
ample, in Model WSWa we obtain the energy per SN in the range 

0.5 < Es^ < 1.5 X 10^^ erg, 

with the average of 0.9 x 10^^ erg. 

The SN sites are randomly distributed in the horizontal coor- 
dinates {x,y). Their vertical positions are drawn from the Gaus- 
sian distributions in z with the scale heights of hu = 0.09 kpc for 
Type II and hi — 0.325 kpc for Type I SNe. Thus, Eq. (1) contains 
the mass source of AMq per SN, 

whereas Eqs. (2) and (3) include kinetic and thermal energy sources 
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of similar strength adding up to ^sn per SN: 

<7SN ^ |SSN + ^) • 

The only other constraints applied when choosing SN sites are to 
reject a site if an SN explosion would result in a local temperature 
above 10^° K or if the local gas number density exceeds 2 cm~^. 
The latter requirement ensures that the thermal energy injected is 
not lost to radiative cooling before it can be converted into kinetic 
energy in the ambient gas. More elaborate prescriptions can be sug- 
gested to select SN sites (Korpi et al. 1999; de Avillez 2000; Joung 
& Mac Low 2006; Gressel et al. 2008); we found this unnecessary 
for our present purposes. 

Arguably the most important feature of SN activity, in the 
present context, is the efficiency of evolution of the SN energy 
from thermal to kinetic energy in the ISM, a transfer that occurs 
via the shocked, dense shells of SN remnants. Given the relatively 
low resolution of this model (and most, if not all, other models 
of this kind), it is essential to verify that the dynamics of expand- 
ing SN shells is captured correctly: inaccuracies in the SN rem- 
nant evolution would indicate that our modelling of the thermal and 
kinetic energy processes was unreliable. Therefore, we present in 
Appendix B detailed numerical simulations of the dynamical evo- 
lution of an individual SN remnant at spatial grid resolutions in 
the range A = 1-4 pc. We allow the SN remnant to evolve from 
the Sedov-Taylor stage (at which SN remnants are introduced in 
our simulations) for t ^ 3.5 Myr. The remnant enters the snow- 
plough regime, with a final shell radius exceeding 100 pc, and we 
compare the numerical results with the analytical solution of Cioffi 
et al. (1998). The accuracy of the numerical results depends on the 
ambient gas density no : larger no requires higher resolution to re- 
produce the analytical results. We show that agreement with Cioffi 
et al. (1998) in terms of the shell radius and expansion speed is 
excellent at resolutions A ^ 2pc for no ^ lcm~^, and also 
very good at A = 4pc for no ~ 0.1 and 0.01 cm~^. Compar- 
isons with models of higher resolution (de Avillez & Breitschwerdt 
2004; Joung et al. 2009), in Section 8.3, also indicate that our basic 
A = 4 pc resolution is adequate. 

Since shock waves in the immediate vicinity of an SN site are 
usually stronger than anywhere else in the ISM, these tests also 
confirm that our handling of shock fronts is sufficiently accurate 
and that the shock-capturing diffusivities that we employ do not 
unreasonably affect the shock evolution. 

Our standard resolution is A = 4pc. To be minimally re- 
solved, the initial radius of an SN remnant must span at least two 
grid points. Because the origin is set between grid points, a mini- 
mum radius of 7 pc for the energy injection site is sufficient. The 
size of the energy injection region in our model must be such that 
the gas temperature is above 10^ K and below 10^ K: at both higher 
and lower temperatures, energy losses to radiation are excessive 
and adiabatic expansion cannot be established. Following Joung 
& Mac Low (2006), we adjust the radius of the energy injection 
region to be such that it contains GOMq of gas. For example, in 
model WSWa this results in a mean rsN of 35 pc, with a standard 
deviation of 25 pc and a maximum of 200 pc. The distribution of 
radii appears approximately lognormal, so rsN > 75 pc is very in- 
frequent and the modal value is about 10 pc; this corresponds to 
the middle of the Sedov-Taylor phase of the SN expansion. Un- 
like Joung & Mac Low (2006), we found that mass redistribution 
within the injection site was not necessary. Therefore we do not im- 
pose uniform site density, particularly as it may lead to unexpected 



Table 1. The cooling function of Wolfireetal.( 1995) at T < 10^K,joined 
to that of Sarazin & White (1987) at higher temperatures, with A = for 
T < 10 K. This cooling function is denoted WSW in the text (and in the 
labels of our numerical models). 



Tk [K] Ak [ergg-2 s'^ cm^ K'^k] 



10 


3.70 X 10^^ 


2.12 


141 


9.46 X IQis 


1.00 


313 


1.18 X 10^0 


0.56 


6102 


1.10 X 10^° 


3.21 


10^ 


1.24 X 10^7 


-0.20 


2.88 X 10^ 


2.39 X 10"^^ 


-3.00 


4.73 X 10^ 


4.00 X 10^6 


-0.22 


2.11 X 10^ 


1.53 X 10"^"^ 


-3.00 


3.98 X 10^ 


1.61 X 10^2 


0.33 


2.00 X 10^ 


9.23 X 10^0 


0.50 



Table 2. The cooling function of Rosen et al. (1993), labelled RBN in the 
text (and in the labels of our numerical models), with A = for T < 10 K. 



Tfc[K] 


Ajt [erg g 2 s -"^ cm^ K ] 




10 


9.88 X 10^ 


6.000 


300 


8.36 X 10^^ 


2.000 


2000 


3.80 X 10^'^ 


1.500 


8000 


1.76 X 10^2 


2.867 


10^ 


6.76 X 10^9 


-0.650 


10^ 


8.51 X 10^2 


0.500 



consequences in the presence of magnetic fields in our MHD sim- 
ulations (described elsewhere). 

2.3 Radiative cooling and photoelectric heating 

We consider two different parameterizations of the optically thin 
radiative cooling appearing in Eq. (3), both of the piecewise power- 
law form A = AkT^^ within a number of temperature ranges 
Tk ^ T < Tjt+i, with Tk and Ak given in Tables 1 and 2. Since 
this is just a crude (but convenient) parameterization of numerous 
processes of recombination and ionisation of various species in the 
ISM, there are several approximations designed to describe the va- 
riety of physical conditions in the ISM. Each of the earlier models 
of the SN-driven ISM adopts a specific cooling curve, often with- 
out explaining the reason for the particular choice or assessing its 
consequences. In this paper, we discuss the sensitivity of the results 
to the choice of the cooling function. 

One parameterization of radiative cooling, labelled WSW and 
shown in Table 1, consists of two parts. For T < 10^ K, we use 
the cooling function fitted by Sanchez- Salcedo et al. (2002) to the 
'standard' equilibrium pressure-density relation of Wolfire et al. 
(1995, cf. Fig. 3b therein). For higher temperatures, we adopt the 
cooling function of Sarazin & White (1987). This part of the cool- 
ing function (but extended differently to lower temperatures) was 
used by Slyz et al. (2005) to study star formation in the ISM. The 
WSW cooling function was also used by Gressel et al. (2008). It has 
two thermally unstable ranges: at 313 ^ T < 6102 K, the gas is 
isobarically unstable (/3k < 1); at T > 10^ K, gas is isochorically 
or isentropically unstable (jSk < and /3k < —1.5, respectively). 

Results obtained with the WSW cooling function are com- 
pared with those using the cooling function of Rosen et al. (1993), 
labelled RBN, whose parameters are shown in Table 2. This cool- 
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102 10^ 10* 10^ 10^ 10^ 10^ 
T [K] 

Figure 1. The cooling functions WSW (solid, black) and RBN (red, dash- 
dotted), with parameters given in Tables 1 and 2, respectively. 



ing function has a thermally unstable part only above 10 K. Rosen 
et al. (1993) truncated their cooling function at T = 300 K. Instead 
of abrupt truncation, we have smoothly extended the cooling func- 
tion down to 10 K. This has no palpable physical consequences 
as the radiative cooling time at these low temperatures becomes 
longer (lOMyr) than other time scales in the model, so that adi- 
abatic cooling dominates. The minimum temperature reported in 
the model of Rosen et al. (1993) is about 100 K. Here, with better 
spatial resolution, the lowest temperature is typically below 50 K. 

We took special care to accurately ensure the continuity of the 
cooling functions, as small discontinuities may affect the perfor- 
mance of the code; hence the values of in Table 1 differ slightly 
from those given by Sanchez- Salcedo et al. (2002). The two cool- 
ing functions are shown in Fig. 1. The cooling function used in 
each numerical model is identified with a prefix RBN or WSW in 
the model label (see Table 3). The purpose of Models RBN and 
WSWb is to assess the impact of the choice of the cooling func- 
tion on the results (Section 8.1). Other models employ the WSW 
cooling function. 

We also include photoelectric heating in Eq. (3) via the stellar 
far-ultraviolet (UV) radiation, T, following Wolfire et al. (1995) 
and allowing for its decline away from the Galactic mid-plane with 
a length scale comparable to the scale height of the stellar disc near 
the Sun (cf. Joung & Mac Low 2006): 

V{z) = roexp(-|z|/300pc) , To = 0.0147 erg g"^ s"\ 

This heating mechanism is smoothly suppressed at T > 2 x 10^ K, 
since the photoelectric effect due to UV photon impact on PAHs 
(Polycyclic Aromatic Hydrocarbons) and small dust grains is im- 
peded at high temperatures (cf. Wolfire et al. 1995). 

2.4 Numerical methods 

2.4.1 The computational domain 

We model a relatively small region within the galactic disc and 
lower halo with parameters typical of the solar neighbourhood. 
Using a three-dimensional Cartesian grid, our results have been 
obtained for a region 1.024 x 1.024 x 2.24 kpc^ in size, with 
1.024 kpc in the radial and azimuthal directions and 1.12 kpc ver- 
tically on either side of the galactic mid-plane. Assuming that the 
correlation length of the interstellar turbulence is /o — 0.1 kpc (see 
Section 6), the computational domain encompasses about 2,000 tur- 
bulent cells, so the statistical properties of the ISM can be reliably 
captured. We are confident that our computational domain is suf- 
ficiently broad to accommodate comfortably even the largest SN 



remnants at large heights, so as to exclude any self-interaction of 
expanding remnants through the periodic boundaries. 

Vertically, our reference model accommodates ten scale 
heights of the cold Hi gas, two scale heights of diffuse Hi (the 
Lockman layer), and one scale height of ionised hydrogen (the 
Reynolds layer). The vertical size of the domain in the reference 
model is insufficient to include the scale height of the hot gas, 
and it would be preferable to consider a computational box of a 
larger vertical size, 2Lz . Indeed, some similar ISM models use a 
vertically elongated computational box with the horizontal size of 
1 kpc X 1 kpc but the top and bottom boundaries at Lz = 10 kpc 
(e.g., de Avillez & Breitschwerdt 2007, and references therein). 
However, the horizontal size of the domain in a taller box may 
need to be increased to keep its aspect ratio of order unity, so as to 
avoid introducing other unphysical behaviour at > 

This constraint arises mainly from the periodic (or sliding pe- 
riodic) boundary conditions in the horizontal planes as they pre- 
clude divergent flows at scales comparable to L±. However, the 
scale of the gas flow unavoidably increases with \z\ because of 
the density stratification. The steady- state continuity equation for 
a gas stratified in 2;, V • i6 = —Uz'd In p/dz, leads to the following 
estimate of the horizontal perturbation velocity arising due to the 
stratification: 



(8) 



where H is the density scale height, dhip/dz —H~^, 
and l± is the horizontal scale of the flow, introduced via 
\dux/dx\, \ duy/dy\ u±/l±. Here we have neglected the ver- 
tical variation of Uz, so that V • u ^ dux/dx + duy/dy: this is 
justified for the hot and warm gas, since their vertical velocities 
vary weakly with z at \z\ > 0.3 kpc (see Fig. 12). Assuming for 
the sake of simplicity that u± is a constant, in Eq. (8), where l±o 
is the horizontal correlation length of u± ai z = 0, we obtain the 
following estimate of the horizontal correlation length 3.t\z\=Lz, 
the top of the domain: 



t± \\z\=L, 



^ + u±t^lo{l + Lz/H) 



where the time available for the expansion is taken as t = Lz/uz, 
lo is the horizontal correlation length of at 2; = 0. We find 
lo 0.1 kpc (Table 5) and H 0.5 kpc (Fig. 19), so that the 
correlation scale of the velocity perturbation at the top and bottom 
boundaries of our domain, Lz ^ I kpc, follows as 



3/0 — 0.3 kpc. 



Indeed, we find the correlation scale of the random flow increases 
to 200-300 pc at 2; = 0.8 kpc (Table 5), so that the diameter of the 
correlation cell, 400-600 pc becomes comparable to the horizontal 
size of the domain, L± = 1 kpc. At larger heights, the periodic 
boundary conditions would suppress the horizontal flows, so that 
that the continuity equation could only be satisfied via an unphys- 
ical increase in the vertical velocity with \z\. In addition, the size 
of SN remnants also increases with as the ambient pressure de- 
creases. Thus, the gas velocity field (and other results) obtained in 
a model with periodic boundary conditions in x and y becomes un- 
reliable at heights significantly exceeding the horizontal size of the 
computational domain. 

The lack of a feedback of the halo on the gas dynamics in the 
disc can, potentially, affect our results. However, we believe that 
this is not a serious problem and, anyway, it would not necessarily 
be resolved by using a taller box of a horizontal size of only 1- 
2 kpc. The gas flow from the halo is expected to be in the form of 
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relatively cool, dense clouds, formed at large heights via thermal in- 
stability or accreted from the intergalactic space (e.g. Wakker & van 
Woerden 1997; Putman et al. 2012). A strong direct (as opposed to 
a long-term) effect of this gas on the multi-phase gas structure in 
the disc is questionable, as it provides just a fraction of the disc's 
star formation rate, 0.1-0. 2M0 yr~^ versus O.S-SM© yr~^ (Put- 
man et al. 2012). Anyway, a taller computational domain would 
not help to include the accreted intergalactic gas in simulations of 
this type. In a galactic fountain, gas returns to the disc at a galacto- 
centric distance at least 3 kpc away from where is starts (Bregman 
1980), and this could not be accounted for in models with tall com- 
putational boxes that are only 1-2 kpc big horizontally. 

In light of these concerns, and since it is not yet possible to ex- 
pand our domain significantly in all three dimensions, we prefer to 
restrict ourselves to a box of height Lz ~ 1 kpc, thus retaining an 
aspect ratio of order unity. This choice of a short box requires great 
care in the choice of vertical boundary conditions (which might also 
introduce unphysical behaviour). We discuss our boundary condi- 
tions in detail in appendix C, but briefly note here that we use mod- 
ified open boundary conditions on the velocity at z = =bL^. These 
conditions allow for both inflow and outflow, and so are to some 
extent capable of simulating gas exchange between the disc and the 
halo, driven by processes within the disc. More specifically, matter 
and energy are free to flow out of and into the computational do- 
main across the top and bottom surfaces if the internal dynamics 
so require. (An inflow occurs when pressure beneath the surface is 
lower than at the surface or in the ghost zones). 

2.4.2 Numerical resolution 

For our standard resolution (numerical grid spacing) Ax = Ay — 
= A = 4pc, we use a grid of 256 x 256 x 560 (exclud- 
ing 'ghost' boundary zones). We apply a sixth-order finite differ- 
ence scheme for spatial vector operations and a third-order Runge- 
Kutta scheme for time stepping. We also investigate one model 
at doubled resolution, A = 2pc, labelled WSWah in Table 3; 
the starting state for this model is obtained by remapping a snap- 
shot from the standard-resolution Model WSWa at t = 600 Myr 
(when the system has settled to a statistical steady state) onto a grid 
512 X 512 X 1120 in size. 

Given the statistically homogeneous structure of the ISM in 
the horizontal directions at the scales of interest (neglecting arm- 
interarm variations), we apply periodic boundary conditions in the 
azimuthal {y) direction. Differential rotation is modelled using the 
shearing-sheet approximation with sliding periodic boundary con- 
ditions (Wisdom & Tremaine 1988) in x, the local analogue of 
cylindrical radius. We apply slightly modified open vertical bound- 
ary conditions, described in some detail in Appendix C, to allow 
for the free movement of gas to the halo without preventing in- 
ward flows at the upper and lower boundaries. In the calculations 
reported here, outflow exceeds inflow on average, and there is a 
net loss of mass from our domain, of order 15% of the total mass 
per Gyr. We do not believe that this slow loss of mass significantly 
affects our results 



2.4.3 Transport coefficients 

The spatial and temporal resolutions attainable impose lower limits 
on the kinematic viscosity and thermal diffusivity which are, 
unavoidably, much higher than any realistic values. These limits 
result from the Courant-Friedrichs-Lewy (CFL) condition which 



requires that the numerical time step must be shorter than the cross- 
ing time over the mesh length A for each of the transport processes 
involved. It is desirable to avoid unnecessarily high viscosity and 
thermal diffusivity. The cold and warm phases have relatively small 
perturbation gas speeds (of order lOkms"^), so we prescribe 
and X to be proportional to the local speed of sound, v — viCs/ci 
and X = XiCs/ci. We ensure that the maximum Reynolds and 
Peclet numbers based on the mesh separation A are always close 
to unity throughout the computational domain (see Appendix C): 
ui ^ 4.2 X 10~^ kms~^ kpc, xi ~ 4.1 x 10~^ kms~^ kpc and 
ci = lkms~^. This gives, for example, x — 0.019 kms~^ kpc 
at T = 10^ K and 0.6 km s"^ kpc at T = 10^ K. Thus, transport 
coefficients are larger in the hot gas where typical temperature and 
perturbation velocity are of order 10^ K and lOOkms"^, respec- 
tively. In all models x — O.liy, i.e., the Prandtl number Pr ^ 10. 
The corresponding fluid Reynolds and Peclet numbers, based on 
the correlation scale of the flow, fall in the range 20-40 in the mod- 
els presented here. 

Numerical handling of the strong shocks widespread in the 
ISM needs special care. To ensure that they are always resolved, 
we include shock-capturing diffusion of heat and momentum, with 
the diffusivities Cx C^^^ respectively, defined as 

^ I c^Ax'^ max5 | V • ii| if V • it < 0, 
^ 1 otherwise 

(and similarly for fi,, but with a coefficient c^), where maxs de- 
notes the maximum value occurring at any of the five nearest mesh 
points (in each coordinate). Thus, the shock-capturing diffusivities 
are proportional to the maximum divergence of the velocity in the 
local neighbourhood, and are confined to the regions of convergent 
flow. Here, = c^^ is a dimensionless coefficient which we have 
adjusted empirically to 10. This prescription spreads a shock front 
over sufficiently many (usually, four) grid points. Detailed test sim- 
ulations of an isolated expanding SN remnant in Appendix B con- 
firm that this prescription produces quite accurate results, particu- 
larly those which are relevant to our goals: most importantly, the 
conversion of thermal to kinetic energy in SN remnants. 

With a cooling function susceptible to thermal instability, ther- 
mal diffusivity x has to be large enough as to allow us to resolve its 
most unstable normal modes: 



7 ''"cool V 27r y 

where /3 is the cooling function exponent in the thermally unstable 
range, TcooI is the radiative coohng time and 7 = 5/3 is the adi- 
abatic index. Figure 4 makes it evident that, in our models, TcooI 
typically exceeds 1 Myr in the thermally unstable regime. Further 
details can be found in Appendix D where we demonstrate that, 
with the parameters chosen in our models, thermal instability is 
well resolved by the numerical grid. 

The shock-capturing diffusion broadens the shocks and in- 
creases the spatial spread of density around them. An undesirable 
effect of this is that the gas inside SN remnants cools faster than 
it should, thus reducing the maximum temperature and affecting 
the abundance of the hot phase. Having considered various ap- 
proaches while modelling individual SN remnants in Appendix B, 
we adopt a prescription which is numerically stable, reduces gas 
cooling within SN remnants, and confines extreme cooling to the 
shock fronts. Specifically, we multiply the term (F — pA)T~^ in 
Eq. (3) by 

$ = exp(-C|VCxl'), (10) 
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Table 3. Selected parameters of the numerical models explored in this paper, named in Column (1). Columns (2)-(3) give input parameters: numerical 
resolution A and initial mid-plane gas number density no- The remaining columns give output parameters: (4) time span over which the models have been in 
steady state (in the units ofr = Lx/uo^rms, the typical horizontal crossing time based on the root-mean- square (r.m.s.) random speed uq given in Column (9) 
and Lx ^ 1 kpc); (5) average kinematic viscosity (i/); (6) average sound speed (cg); (7)-(8) average Reynolds numbers defined at the grid spacing, A, and 
based on the correlation scale of the random flow, Iq ~ 100 pc; (9)-(10) r.m.s. perturbation velocity Urms, and r.m.s. random velocity wo,rms; (H) thermal 
energy density eth; (12) kinetic energy density e^in; and (13) volume fractions fv of cold (C), warm (W) and hot (H) gas at |2;| ^ 200 pc. 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


Model 


A 
[pc] 


no 
1 " 


At 
[r] 


W) 

km kpc 

s 


(Cs) 

[kms-1] 


(Re^> 


(Re) 


Urms 

[kms-1] 


'^0,rms 

[kms-1] 


eth 
kpc^ J 




ekin 
kpc*^ 


/v,C:W:H 
[%] 


WSWa 
WSWah 
RBN 
WSWb 


4 
2 
4 
4 


1.8 
1.8 
2.1 
2.1 


3.9 
0.5 
2.7 
4.0 


0.44 
0.77 
0.24 
0.27 


108 
186 
58 
65 


0.88 
0.85 
1.18 
0.97 


22 
43 
30 
24 


76 
103 
37 
45 


26 
34 
18 
20 


30 
19 
25 
29 


13 
10 
9 
13 


2 : 60 : 38 
3:51:46 
3 : 82 : 15 
3 : 70 : 27 



where Cx is the shock diffusivity defined in Eq. (9). Thus, ^ ~ 1 
almost anywhere in the domain but reduces towards zero in strong 
shocks, where |VCxl^ is large. The value of the additional empir- 
ical parameter, C ^ 0.01, was chosen to ensure numerical stabil- 
ity with minimum change to the basic physics. We have verified 
that, acting together with other artificial diffusion terms, this does 
not prevent accurate modelling of individual SN remnants (see Ap- 
pendix B). 



2.4.4 Initial conditions 

We adopt an initial density distribution corresponding to isothermal 
hydrostatic equilibrium in the gravity field of Eq. (7): 



■ po exp 



ai zi 



a2 z 
2ai Z2 



(11) 



Since our present model does not contain magnetic fields or cos- 
mic rays, which provide roughly half of the total pressure in the 
ISM (the remainder coming from thermal and turbulent pressures), 
we expect the gas scale height to be smaller than that observed. 
Given the limited spatial resolution of our simulations, the corre- 
spondingly weakened thermal instability and neglected self-gravity, 
it is not quite clear in advance whether the gas density used in our 
model should include molecular hydrogen or, alternatively, include 
only diffuse gas. 

We used po = 3.5 x IQ-^'^gcm"^ for models RBN and 
WSWb, corresponding to gas number density, no — 2.1 cm~^ at 
the mid-plane. This is the total interstellar gas density, including 
the part confined to molecular clouds. These models, discussed in 
Section 8.2, exhibit unrealistically strong cooling. Therefore, the 
subsequent models WSWa and WSWah have a smaller amount 
of matter in the computational domain (a 17% reduction), with 
Po = 3.0 X 10~^^gcm~^, or no = 1.8 cm~^, accounting only 
for the atomic gas (see also Joung & Mac Low 2006). 

As soon as the simulation starts, density-dependent heating 
and cooling affect the gas temperature, so it is no longer isothermal 
and p{z) given in Eq. (1 1) is not a hydrostatic distribution. To avoid 
unnecessarily long initial transients, we impose a non-uniform ini- 
tial temperature distribution so as to be near static equilibrium: 



T{z) 



To 



+ z' + 



a2 z 
2ai Z2 



(12) 



where To is obtained from 

r(0) = poA(To) ^ 0.0147 ergg-' s-\ 



The value of To therefore depends on po and the choice of the cool- 
ing function. 



2.5 Models explored 

We considered four numerical models, with relevant input parame- 
ters listed in Table 3, along with some output parameters describing 
the results. The models are labelled with prefix RBN or WSW ac- 
cording to the cooling function used. Angular brackets in Table 3 
denote averages over the whole volume, taken from eleven snap- 
shots (10 for WSWah) within the statistical steady state. The time 
span. At, is given in Column 4, normafised by r = Lx/uo,rms, 
where ito,rms is the root-mean-square random velocity and Lx ~ 
1 kpc is the horizontal size of the computational domain (e.g., 
T ~ 38 Myr in Model WSWa). As ly is set proportional to the 
speed of sound Cs, it is variable and the table presents its aver- 
age value (ly) = vi (cs) /ci, where = 0.004 km kpc and 
ci = 1 km in all models. The numerical resolution is adequate 
when the mesh Reynolds number, ReA = uA/v, does not ex- 
ceed a certain value (typically between 1 and 10) anywhere in the 
domain, where A is the grid spacing (4pc for all models, except 
for Model WSWah, where A = 2 pc). Therefore, we ensure that 
^^max A/v < 5, where i^max is the maximum perturbation velocity 
at any time and any grid point. The indicative values in Table 3 are 
averages of the mesh Reynolds number, (Re a) = (^o / Cs) Aci / vi , 
and the Reynolds number, (Re) = {uo/ Cs)loci/vi. The Reynolds 
number based on the correlation scale of the random flow, lo ^ 
100 pc, is thus 25 times larger than ReA in all models explored 
here except for Model WSWah, where the difference is a factor of 
50. 

The quantities shown in Table 3 have been calculated as fol- 
lows. In Column 9, the r.m.s. perturbation velocity i^rms is derived 
from the total perturbation velocity field which excludes only 
the overall galactic rotation U . In Column 10, the r.m.s. random 
velocity '^0,rms is obtained with the mean flows (ti)^, defined in 
Eq. (6), deducted from u. In Columns 11 and 12, eth = (pe) and 
ekin = {\p^^) ^re the average thermal and kinetic energy den- 
sities, respectively; the latter includes the perturbation velocity u 
and both are normalised to the SN energy E^^. The values of the 
volume fractions of the cold, warm and hot phases (defined in Sec- 
tion 4) near the mid-plane are given in Column 13. 

The reference model, WSWa, uses the WSW cooling function 
but with lower gas density than WSWb, to exclude molecular hy- 
drogen (see Section 3). Model WSWah, which differs from WSWa 
only in its spatial resolution, is designed to clarify the effects of 
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Figure 2. A 3D rendering of (a) temperature T and (b) density n in 
Model WSWa at t = 551 Myr. Cold, dense gas is mostly restricted to 
near the mid-plane, whereas hot gas extends towards the boundaries. To aid 
visualisation of 3D structure, warm gas (10^ < T < 10^ K) in the panel 
(a) and diffuse (n < 10 ~^ cm~^) in panel (b) have high transparency. 
Thus the extreme temperatures or dense structures are emphasised. 



resolution on the results. We also analyze two models which dif- 
fer only in the cooling function, RBN and WSWb, to assess the 
sensitivity of the results to this choice. 



3 THE REFERENCE MODEL 

Model WSWa is taken as a reference model; it has rotation corre- 
sponding to a flat rotation curve with the Solar angular velocity, and 
gas density reduced to exclude that part which would have entered 
molecular clouds. Results for this model were obtained by the con- 
tinuation of the Model WSWb, in which the mass from molecular 
hydrogen had been included: at t ^ 400 Myr, the mass of gas in 
the domain was changed to that of Model WSWa by reducing gas 




0.0 0.1 0.2 0.5 0.4- 0.5 0.5 

t [Gyr] 



Figure 3. Horizontal (xy) averages of (a) the vertical velocity, (b) temper- 
ature and (c) gas density as functions of time for Model WSWa (Model 
WSWb up to 0.4 Gyr). 



density by 15% at every mesh point. The effect of this change of 
the total mass is discussed in Section 8.2. 

Figure 2 shows typical temperature and density distributions 
in this model at t = 551 Myr (i.e., 151 Myr after the restart from 
Model WSWb with reduced density). Supernova remnants appear 
as irregularly shaped regions of hot, dilute gas. A hot bubble break- 
ing through the cold gas layer extends from the mid-plane towards 
the lower boundary, visible as a vertically stretched region in the 
temperature snapshot near the (x, 2;) -face. Another, smaller one can 
be seen below the mid-plane near the (y, ^)-face. Cold, dense struc- 
tures are restricted to the mid-plane and occupy a small part of the 
volume. Very hot and cold regions exist in close proximity. 

Horizontally averaged quantities are shown in Fig. 3 as func- 
tions of z and time for Model WSWb at t < 400 Myr, and 
WSWa at later times, showing the effect of reducing the total mass 
of gas at the transition time. Average quantities may have limited 
physical significance because the multi-phase gas has an extremely 
wide range of velocities, temperatures and densities. For example, 
panel (b) shows that the average temperature near the mid-plane, 
1 2; I < 0.35 pc, is, perhaps unexpectedly, generally higher than that 
at the larger heights. This is due to Type II SN remnants, which 
contain very hot gas with T > 10^ K and are concentrated near 
the mid-plane; even though their total volume is small, they signif- 
icantly affect the average temperature. 

Nevertheless these help to illustrate some global properties of 
the multi-phase structure. Before the system settles into a quasi- 
stationary state at about t = 250 Myr, it undergoes a few large- 
scale transient oscillations involving quasi-periodic vertical mo- 
tions. The period of approximately 100 Myr, is consistent with 
the breathing modes identified by Walters & Cox (2001) and at- 
tributed to oscillations in the gravity field. Gas falling from high 
altitude overshoots the midplane and thus oscillates around it. Tur- 
bulent and molecular viscosities dampen these modes. At later 
times, a systematic outflow develops with an average speed of 
about 100 km ; we note that the vertical velocity increases very 
rapidly near the mid-plane and varies much less at larger heights. 
The result of the reduction of gas density ait ^ 400 Myr is clearly 
visible, as it leads to higher mean temperatures and a stronger and 
more regular outflow, together with a less pronounced and more 
disturbed layer of cold gas. 
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Figure 4. The joint probability density of the gas number density and tem- 
perature, shown for the whole computational domain using 1 1 snapshots of 
Model WSWa in a statistically steady state for 634 ^ t ^ 644 Myr. Con- 
tours of constant cooling time TcooI =10^, 10^ and 10^ yr, are shown to 
clarify the importance of radiative cooling in the model. 



4 THE MULTI-PHASE STRUCTURE 

All models discussed here have a well-developed multi-phase struc- 
ture apparently similar to that observed in the ISM. Since the 
ISM phases are not genuine, thermodynamically distinct phases 
(e.g. Vazquez- Semadeni 2012), their definition is tentative, with 
the typical temperatures of the cold, warm and hot phases usu- 
ally adopted as T :^ 10^ K, 10^-10^ K and 10^ K, respectively. 
However, the borderline temperatures (and even the number of dis- 
tinct phases) can be model-dependent, and they are preferably de- 
termined by considering the results, rather than a priori. Inspec- 
tion of the probability distribution of gas number density and tem- 
perature, displayed in Fig. 4, reveals three distinct concentrations 
at (T[K],n[cm-3]) = (10^10), (10^ 10"^) and (10^ IQ-^). 
Thus, we can confirm that the gas structure in this model can be rea- 
sonably well described in terms of three distinct phases. Moreover, 
we can identify the boundaries between them as the temperatures 
corresponding to the minima of the joint probability distribution at 
about 500 K and 5 x 10^ K. 

The curves of constant cooling time, also shown in Fig. 4, sug- 
gest that the distinction between the warm and hot gas is due to the 
maximum of the cooling rate near T — 10^ K (see also Fig. 1), 
whereas the cold, dense gas, mainly formed by compression (see 
below), closely follows the curve TcooI ~ 10^ yr. 

In Fig. 5, we show the probability distributions of gas number 
density, random velocity, Mach number, thermal and total pressures 
within each phase in Model WSWa. The overlap in the gas density 
distributions (Fig. 5a) is small (at the probability densities of order 
V — 0.1). The ratios of the probability densities near the maximum 
for each phase (mode) are about 100; the modal densities, n ^ 
10~^, 10~^ and 10 cm~^, thus typify the hot, warm and cold gas 
respectively. 

The velocity probability distributions in Fig 5b reveal a clear 
connection between the magnitude of the random velocity of gas 
and its temperature: the r.m.s. velocity in each phase scales with 
its speed of sound. This is confirmed by the Mach number distri- 
butions in Fig. 5c: both warm and hot phases are transonic with 
respect to their sound speeds. The cold gas is mostly supersonic, 
having speeds typically under lOkms"^. The double peak in the 
probability density for the cold gas velocity (Fig 5b) (and the corre- 
sponding extension of the Mach number distribution to > 1) is 
a robust feature, not dependent on the temperature boundary of the 
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Figure 5. The probability distributions of (a) density, (b) random velocity, 
UQ (c) Mach number of random motions uq (defined with respect to the 
local speed of sound), (d) thermal pressure, and (e) total pressure, for each 
phase of Model WSWa, using 1 1 snapshots spanning t = 634 to 644 Myr 
and presented for each phase: cold T < 500 K (black, solid line), warm 
500 ^ T < 5 X 10^ K (blue, dashed), and hot T ^ 5 x 10^ K (red, 
dash-dotted). 



cold gas. This plausibly includes ballistic gas motion in the shells 
of SN remnants, as well as bulk motions of cold clouds at subsonic 
or transonic speed with respect to the ambient warm gas. 

Probability densities of thermal pressure, shown in Fig. 5d, 
are notable for the relatively narrow spread: one order of mag- 
nitude, compared to a spread of six orders of magnitude in gas 
density. Moreover, the three phases have overlapping thermal pres- 
sure distributions, suggesting that the system is in a statistical ther- 
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Figure 6. The probability distributions of gas desnity in model WSWa for 
the cold (black, solid), warm and hot gas. The warm/hot gas have been 
divided into regions \z\ ^ 200 pc (purple/red, dashed/dash-triple-dotted) 
and \z\ > 200 pc (light blue/orange, dash-dotted/long-dashed). The best- 
fit lognormal distributions are shown dotted in matching colour. 



Table 4. Parameters of the lognormal fits to the distribution of gas number 
density n in various phases, where /in and Sn are defined in Eq. (13). 



Phase 


/in [In cm 2] 


Sn [In cm 2] 


cold 


2.02 


0.92 


warm(|z| ^ 200 pc) 


-1.64 


1.47 


warm ( 2; > 200 pc) 


-3.29 


1.47 


warm (total) 


-3.03 


1.47 


hot(|z| ^ 200 pc) 


-5.78 


1.20 


hot(|z| > 200 pc) 


-6.96 


0.77 



mal pressure balance. However, thermal pressure is not the only 
part of the total pressure in the gas, which here includes the tur- 
bulent pressure ^p\u — {u)^]"^, where (u)^, defined in Eq. (6), 
is the mean fluctuation velocity. As shown in Fig. 13, total ki- 
netic energy within the computational domain, associated with ran- 
dom flows, is about a third of the thermal pressure. Correspond- 
ingly, the total pressure distributions in Fig. 5e peaks at about 
4 X 10~^^ dyncm"^ (or ergcm"^), for both the warm and hot 
gas. The cold gas appears somewhat overpressured, with the modal 
pressure at 2 x 10~^^ dyn cm~^, and with some regions under pres- 
sure as high as 10~^^ dyn cm~^. It becomes apparent (see the dis- 
cussion of Fig 7, below) that this is due to both compression by 
transonic random flows and the vertical pressure gradient. All the 
cold gas occupies the higher pressure mid-plane, while the warm 
and hot gas distributions mainly include lower pressure regions 
away from the disc. 

Cold, dense clouds are formed through radiative cooling fa- 
cilitated by compression, which has more importance than in the 
other, hotter phases. The compression is, however, truncated at the 
grid scale of 4 pc, preventing the emergence of higher densities in 
excess of about 10^ cm~^. 

The probability distributions of gas density in Fig. 5a can 
be reasonably approximated by the lognormal distributions, of the 
form 

V{n) = A(/in,Sn) ^ ^^exp f- ^^^^~/-^' ) ■ (13) 

The quality of the fits is illustrated in Fig. 6, using 500 data bins in 
the range 10"'^"^ < n < 10^"^ cm~^; the best-fit parameters are 
given in Table 4. Note that, in making these fits, we have subdivided 
the hot and warm gas into that near the mid-plane {\z\ ^ 200 pc) 
and that at greater heights {\z\ > 200 pc); the former is located in 
the SN active region, strongly shocked with a broad range of den- 
sity and pressure fluctuations, whereas the latter is predominantly 




log P [erg cm"^] 

Figure 7. ProbabiHty distributions for (a) thermal pressure p and (b) total 
pressure P in Model WSWa, for different gas phases: cold (black, solid); 
warm (blue, dashed) and hot (red, dash-triple-dotted) at ^ 200 pc; 
warm (light blue, dotted) and hot (orange, dotted) at > 200 pc. 



the more diffuse and homogeneous gas in the halo. As can be seen 
in Fig. 6, the shape of the probability distribution of the warm gas 
(rather than the position of its maximum) does not var much with 
\z\. Table 4, thus shows the parameters for the warm in the whole 
volume. The lognormal fits satisfy the Kolmogorov-Smirnov test at 
or above the 95% level of significance. For the hot gas fit the KS 
test fails for the total volume. So only the fits for the hot gas split 
by height are included in Table 4. 

The probability densities of thermal and total pressure, dis- 
played in Fig. 7, show that although the thermal pressure of the 
cold gas near the mid-plane is lower than in the other phases the 
total pressures are much closer to balance. The broad probability 
distribution of the cold gas density is consistent with multiple com- 
pressions in shocks. The hot and warm gas pressure distributions 
are also approximately lognormal. The gas at > 200 pc (dotted 
lines) appears to be in both thermal and total pressure balance. 

In summary, we conclude that the system is close to the state 
of statistical pressure equilibrium: the total pressure has similar val- 
ues and similar probability distributions in each phase. Joung et al. 
(2009) also conclude from their simulations that the gas is in both 
thermal and total pressure balance. This could be expected, since 
the only significant deviation from the statistical dynamic equilib- 
rium of the system is the vertical outflow of the hot gas and en- 
trained warm clouds (see Section 7). 

5 THE FILLING FACTOR AND FRACTIONAL VOLUME 
5.1 Filling factors: basic ideas 

ThQ fractional volume of the ISM occupied by the phase i is given 
by 



where Vi is the volume occupied by gas in the temperature range 
defining phase i and V is the total volume. How the gas is dis- 
tributed within a particular phase is described by the phase filling 
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factor 



(15) 



where the over bar denotes a phase average, i.e., an average only 
taken over the volume occupied by the phase i. describes 
whether the gas density of a phase is homogeneous = 1) or 
clumpy < 1). Both of these quantities are clearly important 
parameters of the ISM, allowing one to characterise, as a function 
of position, both the relative distribution of the phases and their in- 
ternal structure. As discussed below, the phase filling factor is also 
directly related to the idea of an ensemble average, an important 
concept in the theory of random functions and so 0^ provides a 
useful connection between turbulence theory and the astrophysics 
of the ISM. Both fv,i and 0^ are easy to calculate in a simulated 
ISM by simply counting mesh-points. 

In the real ISM, however, neither fv nor (f)i can be di- 
rectly measured. Instead the volume filling factor can be derived 
(Reynolds 1977; Kulkami & Heiles 1988; Reynolds 1991), 



(16) 



for a given phase z, where the angular brackets denote a volume 
average, i.e., taken over the total volume. ^ 

Most work in this area to-date has concentrated on the diffuse 
ionised gas (or warm ionised medium) since the emission measure 
of the free electrons EM oc rig and the dispersion measure of pul- 
sars DM oc rie, allowing $ to be estimated along many lines-of- 
sight (e.g. Reynolds 1977; Kulkarni & Heiles 1988; Reynolds 1991; 
Berkhuijsen et al. 2006; Hill et al. 2008; Gaensler et al. 2008). It is 
useful to generalise the tools derived to interpret the properties of 
a single ISM phase for the case of the multiphase ISM, as this can 
help to avoid potential pitfalls when combining data from different 
sources with similar- sounding names (filling factor, filling fraction, 
fractional volume, etc.) but subtly different meanings. In particu- 
lar, only under the very specific conditions explained below, do the 
volume fiUing factors of the different phases of the ISM sum to 
unity. 

In terms of the volume Vi occupied by phase z. 



(17) 



whilst 



- 1 f 

(^0 = ^ / j riidV, (18) 



the final equality holding because = outside the volume Vi 
by definition. Since the two types of averages differ only in the vol- 
ume over which they are averaged, they are related by the fractional 
volume: 



and 



{rii) = —Ui = fv,ini 



/ 2\ Vi 2 p 2 



(19) 



(20) 



As with the density filling factors introduced here, filling factors of tem- 
perature and other variables can be defined similarly to Eqs. (15) and (16), 



Consequently, the volume filling factor $n,i and the phase filling 
factor (f)n^i are similarly related: 



rii 



fv,i^i- 



(21) 



Thus the parameters of most interest, fv,i and (l)n,i, charac- 
terizing the fractional volume and the degree of homogeneity of a 
phase respectively, are related to the observable quantity by 
Eq. (21). This relation is only straightforward when the ISM phase 
can be assumed to be homogeneous or if one has additional sta- 
tistical knowledge, such as the probability density function, of the 
phase. In the next sub-section we use two simple examples to illus- 
trate how the ideas developed here can be appHed to the real ISM; 
we then use them to develop a new interpretation of existing obser- 
vational data and finally discuss how the properties of our simulated 
ISM compare to observations. But first a brief note about different 
methods of averaging is necessary. 

5.1.1 Averaging methods for observations and theory 

An important feature of the definition of the volume filling factor 
given by Eq. (15), is that the averaging involved is inconsistent with 
that used in theory of random functions. In the latter, the calculation 
of volume (or time) averages is usually complicated or impossi- 
ble and, instead, ensemble averages (i.e. averages over the relevant 
probability distribution functions) are used; the ergodicity of the 
random functions is relied upon to ensure that the two averages are 
identical to each other (Section 3.3 in Monin & Yaglom 2007; Ten- 
nekes & Lumley 1972). But the volume filling factors are not 
compatible with such a comparison, as they are based on averag- 
ing over the total volume, despite the fact that each phase occupies 
only a fraction of it. In contrast, the phase averaging used to derive 
01 is performed only over the volume of each phase, and so should 
correspond better to results from the theory of random functions. 

5.2 Filling factors: applications 

5.2.7 Assumption of homogeneous phases 

The simplest way to interpret an observation of the volume filling 
factor $i is to assume that each ISM phase has a constant density. 
Consider Eqs. (14), (15) and (16) for an idealised two-phase sys- 
tem, where each phase is homogeneous. (These arguments can eas- 
ily be generalised to an arbitrary number of homogeneous phases.) 
For example a set of discrete clouds of one phase, of constant den- 
sity and temperature, embedded within the other phase, with differ- 
ent (but also constant) density and temperature. Let one phase have 
(constant) gas number density A^i and occupy volume Vi , and the 
other N2 and V2, respectively. The total volume of the system is 
y = Vi + y2. 

The volume- averaged density of each phase, as required for 
Eq. (16), is given by 

(n,) = ^ = fy^^N.. (22) 

where i = 1,2. Similarly, the volume average of the squared den- 
sity is 



V 



= fv,iNl 



The fractional volume of each phase can then be written as 



for example (f)T,i = Ti /Tf, etc. 



fv,i = 



jrii)'^ _ (rii) 
{nl) ~ Ni 



(23) 



(24) 
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with fv,i -\- fv,2 = 1, and $i + $2 = 1. The volume- averaged 
quantities satisfy (n) = (ni) + (722) = fv,iNi + /v,2A^2 and 
(n^) = (n?) + (nl) = /v,iA^i 2 iV|, with the density variance 
cr" = (n^) - (n)' = fv,ifvANi - N2Y . 

In contrast, the phase-averaged density of each phase, used to 
calculate the phase filling factor Eq. (15), is simply 7i~ = Ni, and 
the phase average of the squared density is nf = Nf, so that the 
phase filling factor is = 1, as must be the case for a homoge- 
neous phase. 

Thus for homogeneous phases, the volume filling factor and 
the fractional volume of each phase are identical to each other, 
= fv,i, and both sum to unity when considering all phases; in 
contrast, the phase- averaged filling factor is unity for each phase, 
0^ = 1. If a given phase occupies the whole volume (i.e., we have 
a single-phase medium), then all three quantities are simply unity: 
(j)i=^i= fv,i = 1. 

Whilst an assumption of homogeneous phases may be justified 
for some ISM phases, perhaps in specific regions of the galactic 
disc, in the case of the simulated ISM discussed in this paper such 
an assumption would lead to significant underestimates of fv,i for 
all phases, by a factor of 2 for the cold and hot gas and by an order 
of magnitude for the warm gas. 



5.2.2 Assumption of lognormal phases 

For the more realistic case of an inhomogeneous ISM, where each 
phase consists of gas with a range of densities, the interpretation of 
$i requires additional knowledge about the statistical properties of 
a phase. 

For electrons in the diffuse ionised gas Reynolds (1977) de- 
rived the correction factor al/nl^ where ric is the average density 
of electron clouds and cr^ the density variance within clouds, to 
allow for dumpiness in the electron distribution when calculating 
the fraction of the total path length occupied by the clouds. More 
generally, the probability distribution function of the gas in a phase 
allows 0i to be calculated directly, as we now illustrate for the case 
of the lognormal PDFs identified in Section 4. 

For a lognormal distribution V(ni) ^ A(/ii, s^), Eq. (13), the 
mean and mean-square densities are given by the following phase 
('ensemble') averages: 

ni = e^^+*^/^, ai = {ni — nlY = nf {e^^ — 1^ , (25) 
where erf is the density variance around the mean ru, so that 

—2 —2 

0n,i = = = 2 1'— 2 = exp(-s- ). (26) 
(^i -rrii 

So the phase filling factor 0^,^ = 1 only for a homogeneous density 
distribution, ai = (or equivalently, Si — 0). This makes it clear 
that this filling factor, defined in terms of the phase average, is quite 
distinct from the fractional volume, fv,i, but rather quantifies the 
degree of homogeneity of the gas distribution within a given phase. 
Both describe distinct characteristics of the multi-phase ISM, and, 
if properly interpreted, can yield rich information about the struc- 
ture of the ISM. 

In the case of the simulated ISM, using the lognormal descrip- 
tion of the phases given in Table 4 gives a reasonable agreement 
between the actual and estimated fv,i and (t)i for all phases, with 
the biggest discrepancy being an underestimate of fv^warm ~ 0.4 
against a true value of /A/,warm ~ 0.6. 



5.2.3 Application to observations 

Observations can be used to estimate the volume-averaged filling 
factor defined in Eq. (16), for a given ISM phase. On its own, 
this quantity is of limited value in understanding how the phases 
of the ISM are distributed: of more use are the fractional volume 
occupied by the phase fv,i, defined in Eq. (14), and its degree of 
homogeneity which is quantified by defined by Eq. (15). Know- 
ing $i and fv^i follows via Eq. (21): 

fv,i = (27) 

This formula is exact, but its applicability in practise is limited if 
(pi is unknown. However 0^ can be deduced from the probability 
distribution of rii : for example if the the density probability distri- 
bution of the phase can be approximated by the lognormal, as is 
expected for a turbulent compressible gas (Vazquez-Semadeni & 
Garcia 2001; Elmegreen & Scalo 2004), then can be estimated 
from Eq. (26). 

To illustrate how these quantities may be related, let us con- 
sider some observations reported for the diffuse ionised gas (the 
general approach suggested can be applied to any observable or 
computed quantity). Berkhuijsen et al. (2006) and Berkhuijsen & 
Miiller (2008) estimated $dig for the diffuse ionised gas (DIG) 
in the Milky Way using dispersion measures of pulsars and emis- 
sion measure maps. In particular, Berkhuijsen et al. (2006) obtain 
^DiG — 0.24 towards \z\ — Ikpc, and Berkhuijsen & Miiller 
(2008) find the smaller value $dig — 0.08 for a selection of pul- 
sars that are closer to the Sun than the sample of Berkhuijsen et al. 
(2006). On the other hand, Berkhuijsen & Fletcher (2008, 2012) 
used the same data for pulsars with known distances to derive PDFs 
of the distribution of average DIG cloud densities which are well 
described by a lognormal distribution; the fitted lognormals have 
SDiG ^ 0.32 (Table 1 in Berkhuijsen & Fletcher 2012). Using 
Eqs. (26) and (27), this implies that the fractional volume of DIG 
with allowance for its inhomogeneity is about 

/V,DIG ^ 0.1-0.3. 

In other words the combination of $dig and sdig from these 
results imply that the DIG is approximately homogeneous. This 
value of /v,DiG is in good agreement with the earher estimates 
of Reynolds (1977) and Reynolds (1991) who obtained /f,dig ^ 
0.1-0.2 and close to that of Hill et al. (2008) who obtained 
/f,dig ~ 0.25 for a vertically stratified ISM, by comparing ob- 
served emission and dispersion measures to simulations of isother- 
mal MHD turbulence. 

Volume density PDFs derived from observations are still rare. 
However, PDFs of the column density (and similar observables 
such as emission measure and dispersion measure) are more eas- 
ily derived. The applicability of the method outhned in this Sec- 
tion, of deriving the fractional volume occupied by different ISM 
phases from the (observable) volume filling factor and the PDF of 
the density distribution, would improve as the relation between the 
statistical parameters of volume and column density distributions 
becomes better understood. 



5.2.4 Simulation results 

The filling factors and fractional volumes from Equations (14), 
(16) and (15) have been computed for the phases identified in Sec- 
tion 4 for the reference model WSWa and presented in Fig. 8. Vol- 
umes are considered as discrete horizontal slices. To isolate the z- 
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Figure 8. Vertical profiles of (a) the phase-averaged density filling factors 
4>i = /"iT^i of the gas phases identified in the text: cold (black, solid line, 
T < 500 K); warm (blue, dashed, 5 x 10^ ^ T < 5 x 10^ K) and hot 
(red, dash-dotted, T ^ 5 x 10^ K); and (b) the volume-averaged density 
filling factors = {rii)'^ / {nf), and (c) the fractional volumes fv,i with 
the same line style for each phase. The various filling factors are defined and 
discussed in Section 5. These results are from 21 snapshots in the interval 
636 ^ t ^ 646 Myr for Model WSWa. 



dependence we averaged over slices of single cell thickness (4 pc- 
thick). 

The hot gas (Fig. 8c) accounts for about 70% of the volume 
at l^l 1 kpc and about 40% near the mid-plane. The local maxi- 
mum of the fractional volume of the hot gas at < 200 pc is due 
to the highest concentration of SN remnants there, filled with the 
very hot gas. Regarding its contribution to integrated gas parame- 
ters, it should perhaps be considered as a separate phase. 

At ^ 0.7 kpc the warm gas accounts for over 50% of the 
volume. The cold gas occupies a negligible volume, even in the 
mid-plane where it is concentrated. It is, however, quite homoge- 
neous at low \z\ compared to the warm and hot phases, which only 
become relatively homogeneous at \z\ > 0.3 kpc (Fig. 8a). 




200 300 
|/| [pc] 



Figure 9. The second-order structure functions calculated using Eq. (28), 
for the layer —10<z< 10 pc, of the velocity components Ux (black, 
solid line), Uy (blue, dashed) and Uz (red, dash-dot). The offset / is confined 
to the (x, y) -plane only. 



6 THE CORRELATION SCALE OF THE RANDOM 
FLOWS 

We have estimated the correlation length of the random velocity u 
at a single time step of the model WSWa, by calculating the second- 
order structure functions V{1) of the velocity components Ux, Uy 
and Uz , where 



V{l) = ([u{x + l)-u{x)f), 



(28) 



with X the position in the (x, 2/) -plane and I a horizontal offset. 
We did not include offsets in the z-direction and aggregated the 
squared differences by |Z| only. Since the flow is expected to be sta- 
tistically homogeneous horizontally, while the correlation length is 
expected to vary with z. A future paper will analyse in more detail 
the three-dimensional properties of the random flows, including its 
anisotropy and dependence on height. We measured V{1) for five 
different heights, z = 0, 100, -100, 200 pc and 800 pc, averaging 
over six adjacent slices in the (x,^) -plane at each position, cor- 
responding to a layer thickness of 20 pc. The averaging took ad- 
vantage of the periodic boundaries in x and y; for simplicity we 
chose a simulation snapshot at a time for which the offset in the 
^-boundary, due to the shearing boundary condition, was zero. The 
structure function for the mid-plane (— 10 < z < 10 pc) is shown 
in Fig. 9. 

The correlation scale can be estimated from the form of the 
structure function since velocities are uncorrected if I exceeds the 
correlation length /o, so that V becomes independent of /, V{1) ^ 
2i^rms for I ^ Iq. Precisely which value of P(/) should be chosen 
to estimate lo in a finite domain is not always clear; for example, 
the structure function of Uy in Fig. 9 allows one to make a case 
for either the value at which P(/) is maximum or the value at the 
greatest /. Alternatively, and more conveniently, one can estimate 
lo via the autocorrelation function C(/), related to V{1) by 



C{1) = 1 



(29) 



In terms of the autocorrelation function, the correlation scale , 
defined as 



-r 

Jo 



C{l)dl, 



(30) 



and this provides a more robust method of deriving lo in a finite 
domain. Of course, the domain must be large enough to make C{1) 



14 Gent et al. 



Table 5. The correlation scale /q and rms velocity Wrms at various distances 
from the mid-plane. 



z 




[kms 






/o [pc] 




Ux 


Uy 


Uz 
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Uz 
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37 


99 
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50 
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95 


87 


171 


200 


27 


20 


63 


119 


105 


186 


800 


51 


21 


107 


320 


158 


277 



negligible at scales of the order of the domain size; this is a non- 
trivial requirement, since even an exponentially weak tail can make 
a finite contribution to /o • In our estimates we are, of course, lim- 
ited to the range of within our computational domain, so that 
the upper limit in the integral of Eq. (30) is equal to Lx — Ly, the 
horizontal box size. 

Figure 10 shows C(V) for five different heights in the disc, 
where i^rms was taken to correspond to the absolute maximum of 
the structure function, i^^ms = \ niax(P), from Eq. (29) at each 
height. 

The autocorrelation function of the vertical velocity varies 
with z more strongly than, and differently from, the autocorrelation 
functions of the horizontal velocity components; it broadens as |z| 
increases, meaning that the vertical velocity is correlated over pro- 
gressively greater horizontal distances. Already at |z| ^ 200 pc, 
Uz is coherent across a significant horizontal cross-section of the 
domain, and at 800 pc so is Ux- An obvious explanation for 

this behaviour is the expansion of the hot gas streaming away from 
the mid-plane, which thus occupies a progressively larger part of 
the volume as it flows towards the halo. 

Table 5 shows the rms velocities derived from the structure 
functions for each component of the velocity at each height, and 
the correlation lengths obtained from the autocorrelation functions. 
Note that these are obtained without separation into phases. The un- 
certainties in i^rms due to the choices of local maxima in P(/) are 
less than 2 km s~^. However, these can produce quite large system- 
atic uncertainties in /q, as small changes in iXrms can lead to C(X) 
becoming negative in some range of / (i.e. a weak anti-correlation), 
and this can significantly alter the value of the integral in Eq. (30). 
Such an anti-correlation at moderate values of / is natural for in- 
compressible flows; the choice of i^rms and the estimate of /o are 
thus not straightforward. Other choices of i^rms in Fig. 9 can lead 
to a reduction in /q by as much as 30 pc. Better statistics, derived 
from data cubes for a number of different time-steps, will allow for 
a more thorough exploration of the uncertainties, but we defer this 
analysis to a later paper. 

The rms velocities given in Table 5 are compatible with the 
global values of i^rms and for the reference Model WSWa shown 
in Table 3. The increase in the rms value of Uz with height, from 
about 40kms"^ at 2; = to about 60kms"^ at z = 200 pc, 
reflects the systematic net outflow with a speed increasing with \z\. 
There is also an apparent marginal tendency for the rms values of 
Ux and Uy to decrease with increasing distance from the mid-plane. 

The correlation scale of the random flow is very close to 
100 pc in the mid-plane, and we have adopted this value for /o 
elsewhere in the paper. This estimate is in good agreement with 
the hydrodynamic ISM simulations of Joung & Mac Low (2006), 
who found that most kinetic energy is contained by fluctuations 
with a wavelength (i.e. 2/o in our notation) of 190 pc. In the MHD 
simulations of Korpi et al. (1999), /o for the warm gas was 30 pc 
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Figure 10. Autocorrelation functions for the velocity components Ux 
(black, solid line), Uy (blue, dashed) and Uz (red, dash-dot) for 20 pc thick 
layers centred on four different heights, from top to bottom: — 10 < 2; < 
10 pc, 90 < z < 110 pc, -110 < z < -90 pc, 190 < z < 210 pc and 
790 < z < 810 pc. 
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Figure 11. The perturbation velocity field u in Model WSWa at t = 
550 Myr. The colour bar indicates the magnitude of the velocity field de- 
picted in the volume shading, with rapidly moving regions highlighted with 
shades of red. The low velocity regions, shaded blue, have reduced opac- 
ity to assist visualisation. Arrow length of vectors (a) is proportional to the 
magnitude of u, with red (blue) arrows corresponding to Uz > (uz < 0) 
and independent of the colour bar. Trajectories of fluid elements (b) are also 
shown, indicating the complexity of the flow and its pronounced vortical 
structure. 



at all heights, but that of the hot gas increased from 20 pc in the 
mid-plane to 60pc at = 150 pc. de Avillez & Breitschwerdt 
(2007) found lo = 73 pc on average, with strong fluctuations in 
time. As in Korpi et al. (1999), there is a weak tendency for Iq of 
the horizontal velocity components to increase with 1 2; | in our sim- 
ulations, but this tendency remains tentative, and must be examined 
more carefully to confirm its robustness. 
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Figure 12. Contours of the probabiHty density of the vertical velocity Uz 
as a function of z in Model WSWa from 1 1 snapshots at t = 634-644 Myr. 
The cold (T < 500 K), warm (500 K ^ T < 5 x 10^ K) and hot (T ^ 
5 X 10^ K) are shown in panels (a) to (c), respectively. The horizontal 
averages of the vertical velocity Uz in each case are shown in red, dashed 
in each panel as well as the mid-plane position (black, dotted). 



7 GAS FLOW TO AND FROM THE MID-PLANE 

Figure 1 1 illustrates the 3D structure of the perturbation velocity 
field for the reference Model WSWa. Shades of red show the re- 
gions of high speed, whereas regions moving at speeds below about 
300kms~^ are transparent to aid visualisation. Velocity vectors 
are shown in panel (a) using arrows, with size indicating the speed, 
and colour indicating the sign of the ^-component of the velocity 
(indicating preferential outflow from the mid-plane). Red patches 
are indicative of recent SN explosions, and there is a strongly di- 
vergent flow close to the middle of the xz-face. In addition, stream 
lines in panel (b) display the presence of considerable small scale 
vortical flow near the mid-plane. 

The mean vertical flow is dominated by the high velocity 
hot gas, so it is instructive to consider the velocity structure of 
each phase separately. Figure 12 shows the probability distributions 
V{z, Uz) as functions of Uz in the {z, Uz)-p\^ne from 11 snapshots 
of Model WSWa, separately for the cold (a), warm (b) and hot gas 
(c). The cold gas is mainly restricted io\z\ < 300 pc and its vertical 
velocity varies within ±20 km . As indicated by the red dashed 
curve in Panel (a), on average, the cold gas moves towards the mid- 
plane, presumably after cooling at larger heights. The warm gas 
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Figure 13. Evolution of the volume-averaged thermal energy density 
(black: model WSWb, blue: model WSWa, purple: model WSWah, red: 
model RBN) and kinetic energy density (as above; lower lines) in the sta- 
tistical steady regime, normalised to the SN energy ^^sNkpc"^. Mod- 
els WSWb (black) and RBN (red) essentially differ only in the choice of 
the radiative cooling function. 



is involved in a weak net vertical outflow above \z\ = 100 pc, 
of order ±10 km s~^. This might be an entrained flow within the 
hot gas. However, due to its skewed distribution, the modal flow 
and thus mass transfer is typically towards the mid-plane. The hot 
gas has large net outflow speeds, accelerating to about 100 km 
within 1 2; I lb 200 pc, but with small amounts of inward flowing 
gas at all heights. The mean hot gas outflow speed increases at an 
approximately constant rate to somewhat over lOOkms"^ within 
±100 pc of the mid-plane, and then decreases with further distance 
from the mid-plane, at a rate that gradually decreases with height 
for l^l > 0.5 kpc. This is below the escape velocity in the gravita- 
tional potential adopted. The structure of the velocity field shall be 
investigated further elsewhere. 



8 SENSITIVITY TO MODEL PARAMETERS 
8.1 The cooling function 

We consider two models, RBN and WSWb, with parameters given 
in Table 3, to assess the effects of the specific choice of the cool- 
ing function. Apart from different parameterizations of the radia- 
tive cooling, the two models share identical parameters, except the 
value of To was slightly higher in Model RBN, because of the 
sensitivity of the initial conditions to the cooling function (Sec- 
tion 2.4.4). 

The volume- averaged thermal and kinetic energy densities, 
the latter excluding the imposed shear flow U, are shown in 
Fig. 13 as functions of time. The averages for each are shown in 
Columns (11) and (12), respectively of Table 3, using the appropri- 
ate steady state time intervals given in Column (4). Models reach a 
statistical steady- state, with mild fluctuations around a well defined 
mean value, very soon (within 60Myr of the start of the simula- 
tions). The effect of the cooling function is evident: both the ther- 
mal and kinetic energies in Model RBN are about 60% of those 
in Model WSWb. This is understandable as Model RBN has a 
stronger cooling rate than Model WSWb, only dropping below the 
WSW rate in the range T < 10^ K (see Fig. 1). Interestingly, both 
models are similar in that the thermal energy is about 2.5 times the 
kinetic energy. 

These results are also remarkably consistent with results by 
Balsara et al. (2004, their Fig. 6) and Gressel (2008, Fig. 3.1). Gres- 
sel (2008) applies WSW cooling and has a model very similar to 
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Figure 14. Probability density distributions in the whole computational 
domain, obtained without separation into distinct phases, for (a) gas density, 
(b) temperature and (c) thermal pressure, for Model RBN (blue, dashed) and 
Model WSWb (black, solid), in a statistical steady state, each averaged over 
21 snapshots spanning 20Myr (RBN: 266 to 286 Myr, and WSWb: 305 to 
325 My r) and the total simulation domain \z\ ^ 1.12 kpc. The smaller 
frames to the right display the same information but near the midplane, 
\z\ < 20 pc, only. 



Model WSWa, with half the resolution and \z\ ^ 2 kpc. He re- 
ports average energy densities of 24 and lO^sNkpc"^ (thermal 
and kinetic, respectively) with SN rate = asN, comparable to 30 
and 13 ^sn kpc~^ obtained here for Model WSWa. 

Balsara et al. (2004) simulate an unstratified cubic region 
200 pc in size, driven at SN rates of 8, 12 and 40 times the Galac- 
tic rate, with resolution more than double that of Model WSWa. For 
SN rates 12(7sn and SdsN, they obtain average thermal energy den- 
sities of about 225 and 160 Esn kpc~^, and average kinetic energy 
densities of 95 and 60 ^sn kpc~^, respectively (derived from their 
energy totals divided by the [200 pc]^ volume). 

To allow comparison with our models, where the SNe energy 
injection rate is lasN, if we divide their energy densities by 12 and 
8, respectively, the energy densities would be 19 and 20 ^sn kpc~^ 
(thermal), and 8 and 7.5 ^sn kpc~^ (kinetic). These are slightly 
lower than our results with RBN cooling (25 and 9 ^sn kpc~^), 
but are below those with WSW (30 and 13 ^sn kpc"^ for WSWa, 
as given above). Balsara et al. (2004) used an alternative cooling 
function (Raymond & Smith 1977), so allowing for some additional 
uncertainty over the net radiative energy losses, the results appear 
remarkably consistent. 

While cooling and resolution may marginally affect the mag- 
nitudes, it appears that thermal energy density may consistently be 
expected to be about 2.5 times the kinetic energy density, in these 
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Figure 15. Probability densities for various variables in individual phases, for Model WSWb (left-hand column of panels) and Model RBN (right-hand 
column): (a) and (f) for gas density; (b) and (g) for random velocity uq; (c) and (h) for the Mach number of the random velocity defined with respect to the 
local sound speed; (d) and (i) for thermal pressure; and (e) and (j) for the total pressure. The cold phase spans T < 500 K (black, solid), the warm gas has 
500 < T < 5 X 10^ K (blue, dashed) and the hot gas is at T ^ 5 x 10^ K (red, dash-dotted). Eleven snapshots have been used for averaging, spanning 
t = 200-300 Myr for Model RBN and t = 300-400 Myr for Model WSWb. 



models. It also appears, by comparing the stratified and unstratified 
models, that the ratio of thermal to kinetic energy is not strongly 
dependent on height over the range included in our model. 

The two models are further compared in Fig. 14, where we 
show probability distributions for the gas density, temperature and 
thermal pressure. With both cooling functions, the most probable 
gas number density is around 3 x 10~^ cm~^; the most probable 
temperatures are also similar, at around 3 x 10^ K. With the RBN 
cooling function, the density range extends to smaller densities than 
with WSWb; and yet the temperature range for WSWb extends to 



lower values than for RBN. It is evident that the isobarically un- 
stable part of the WSW cooling function does significantly reduce 
the amount of gas at T = 313-6102 K (the temperature range cor- 
responding to the thermally unstable regime of the WSW cooling), 
and increase the amount of gas below 100 K. However this is not 
associated with higher densities than when using the RBN cooling 
function. This may indicate that multiple compressions, rather than 
thermal instability, dominate the formation of dense clouds. 

The most probable thermal pressure is lower in Model RBN 
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than in WSWb, consistent with the lower thermal energy content of 
the former. 

The probability distributions of various quantities, shown in 
Fig. 15, confirm the clear phase separation in terms of gas density 
and perturbation velocity. Here we used the same borderline tem- 
peratures for individual phases as for Model WSWa (Fig. 5). De- 
spite minor differences between the corresponding panels in Figs. 5 
and 15, the peaks in the gas density probability distributions are 
close to 10\ 3 X 10"^ and 10"^ cm~^ in all models. Given the 
extra cooling of hot gas and reduced cooling of cold gas with the 
RBN cooling function, more of the gas resides in the warm phase 
in Model RBN. The thermal pressure distribution in the hot gas re- 
veals the two 'types' (see the end of section 4), which are mostly 
found within \z\ < 200 pc (high pressure hot gas within SN rem- 
nants) and outside this layer (diffuse, lower pressure hot gas). The 
probability distribution for the Mach number in the warm gas ex- 
tends to higher values with the RBN cooling function, perhaps be- 
cause more shocks reside in the more widespread warm gas, at the 
expense of the cold phase. It is useful to remember that, although 
each distribution is normalised to unit underlying area, the frac- 
tional volume of the warm gas is about a hundred times that of the 
cold phase. 

The probability distributions of density and pressure in with- 
out preliminary separation into phases, presented in Fig. 14 do not 
show clear separations into phases (cf. e.g. Joung & Mac Low 
2006; de Avillez & Breitschwerdt 2004), such that division into 
three phases would arguably only be conventional, if based on these 
alone. The probability distributions near the mid-plane, \ z\ < 20 pc 
Fig. 14, exhibit a marginally better phase separation for the gas den- 
sity (smaller frames in Fig. 14) (see also Korpi et al. 1999; Hill et al. 
2012, their Figs. 1; and 6, respectively). However our analysis in 
terms of phase- wise PDFs confirms that the trimodal structure evi- 
dent in the temperature distribution (Fig. 14b) has a complementary 
structure in the gas density. 

Stratification of the thermal structure is clarified in Fig. 16, 
where we introduce narrower temperature bands specified in Ta- 
ble 6. The fractional volume of gas in each temperature range i at a 
height z is given by 



V^izl ^ N,{z) 
V{z) N{z) ' 



(31) 



similarly to Eq. (14), where Ni{z) is the number of grid points 
in the temperature range T^,min ^ T < Ti,max, with Ti,min and 
7i,max given in Table 6, and N{z) is the total number of grid points 
at that height. 

The fractional volumes in Column (13) of Table 3 show that 
near the mid-plane cold gas forms in similar abundances, inde- 
pendent of the cooling function. However, much less hot gas is 
achieved for Model RBN. Figure 16 also helps show how the ther- 
mal gas structure depends on the cooling function. Model WSWb, 
panel (b), has significantly more very cold gas (T < 50 K) than 
RBN, panel (a), but slightly warmer cold gas (T < 500 K) is more 
abundant in RBN. The warm and hot phases (T > 5 x 10^ K) have 
roughly similar distributions in both models, although Model RBN 
has less of both phases. Apart from relatively minor details, the ef- 
fect of the form of the cooling function thus appears to be straight- 
forward and predictable: stronger cooling means more cold gas and 
vice versa. What is less obvious, however, is that the very hot gas 
is more abundant near ±1 kpc in Model RBN than in WSWb, indi- 
cating that the typical densities must be much lower. This, together 
with the greater abundance of cooler gas near the mid-plane, sug- 
gest that there is less stirring with RBN cooling. 



(a) 
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Figure 16. Vertical profiles of the fractional volumes occupied by the var- 
ious temperature ranges, with the key shown in Table 6. (a) Model RBN, 
using 21 snapshots spanning 266 to 286 Myr. (b) Model WSWb, using 21 
snapshots spanning 305 to 325 Myr. 



Table 6. Key to Figs 16 and 17, defining the gas temperature bands used 
there, and the classification within three phases. 



Temperature band 



Line style Phase 



T < 5 X IQi K cold 

5x101K^T<5x102k cold 

5x102k^T<5x103k warm 

5x103k^T<5x104K warm 

5 X 10^ K ^ T < 5 X 10^ K warm 

5 X 10^ K ^ T < 5 X 10^ K — • • — hot 

T ^ 5 X 10^ K hot 



Altogether, we conclude that the properties of the cold and 
warm phases are not strongly affected by the choice of the cool- 
ing function. The main effect is that the RBN cooling function pro- 
duces less hot gas with significantly lower pressures. This can read- 
ily be understood, as this function provides significantly stronger 
cooling at T > 10^ K. 



8.2 The total gas mass 

Models RBN and WSWb have about 17% more mass of gas than 
the reference Model WSWa, where we have removed that part 
of the gas mass which should be confined to molecular clouds 
unresolved in our simulations (as described in section 3). The 
difference is apparent in comparing Fig. 16b with Fig. 17b (or 
Fig. 17a). Higher gas mass causes the abundance of hot gas to re- 
duce with height, contrary to observations, and to the behaviour of 
Model WSWa. Otherwise, the fractional volumes within ±200 pc 
of the mid-plane appear independent of the gas mass. 
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Figure 17. Vertical profiles of the fractional volumes (Eq. 14) for 
Model WSWah (a), which differs only in its doubled spatial resolution 
from the reference Model WSWa (c) and the fractional mass (Eq. 32) from 
Model WSWa (b). These are calculated for the temperature ranges given, 
along with the figure legend, in Table 6. The former uses 10 snapshots and 
the latter 6 spanning 633 to 638 Myr. 



8.3 Numerical resolution 

Models WSWa and WSWah differ only in their resolution, using 2 
and 4 pc, respectively. Model WSWah is a continuation of the state 
of WSWa after 600 Myr of evolution. 

The most obvious effect of increased resolution is the in- 
crease in the magnitude of the perturbed velocity and temperatures; 
('^rms) — 76kms~^ in Model WSWa increasing to 103kms~^ 
in Model WSWah (Table 3, Column 9) and (cs) from 150 to 
230kms~^ (Column 6). Both (u^^^) and the random velocity 
(i^Crms) increased by a similar factor of about 1.3. However, 
the thermal energy eth is reduced by a factor of 0.6 with the higher 
resolution, while kinetic energy ck remains about the same. This 
suggests that in the higher-resolution model, the higher velocities 
and temperatures are associated with lower gas densities. 

The vertical distribution of the fractional volume in each 
temperature range (defined in Table 6) is shown in Fig. 17 for 
Model WSWah (panel a) for comparison with Model WSWa in 




logn [cm~^] 




-13 -12 

logp [erg cm"^] 



Figure 18. Volume weighted probability distributions of gas number den- 
sity (a), temperature (b) and thermal pressure (c) for models WSWa (black, 
solid) and WSWah (blue, dashed) for the total numerical domain \z\ ^ 
1.12 kpc. 



(c). The fractional mass (b) is calculated similarly to Eq. (31) is 



fM,i{z) 



M{z) ' 



(32) 



where Mi{z) is the mass of gas within temperature range i at a 
given z, and M{z) is the total gas mass at that height. 

Note that the relative abundances of the various phases in these 
models might be affected by the unrealistically high thermal con- 
ductivity adopted. The coldest gas (black, solid), with T < 50 K, 
is largely confined within about 200 pc of the mid-plane. Its frac- 
tional volume (Fig. 17a,c) is small even at the mid-plane, but it pro- 
vides more than half of the gas mass at z = (Fig. 17b). Gas in the 
next temperature range, 50 < T < 500 K (purple, dash-dotted), is 
similarly distributed in z. Models WSWa and WSWah differ only 
in their resolution, using 2 and 4 pc, respectively. Model WSWah 
is a continuation of the state of WSWa after 600 Myr of evolution. 
With higher resolution the volume fraction of the coldest gas is 
significantly enhanced (Fig. 17c compared to a), but it is similarly 
distributed. 

Gas in the range 5 x 10^ < T < 5 x 10^ K (dark blue, 
dashed) has a similar profile to the cold gas for both the fractional 
mass and the fractional volume, and this is insensitive to the model 
resolution. This is identified with the warm phase, but exists in the 
thermally unstable temperature range. It accounts for about 10% by 
volume and 20% by mass of the gas near the mid-plane, which is 
consistent with observational evidence. It is negligible away from 
the supernova active regions. 

The two bands with T > 5 x 10^ K (red, dotted and or- 
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Figure 19. Horizontal averages of gas number density, n(z) (a), and total 
pressure, P(z) (b), for Model WSWa (solid, black), and Model WSWah 
(dashed, blue). Each are time-averaged using 6 and 10 snapshots respec- 
tively, spanning 633 to 638 Myr. The vertical lines indicate standard devia- 
tion within each horizontal slice. The thermal p{z) (dotted) and ram po{z) 
(fine dashed) pressures are also plotted (b). 



ange, dash- 3 dotted) behave similarly to each other (Fig. 17a,c), 
occupying similar fractional volumes for |z| < 0.75 kpc, and 
with fv,i increasing above this height (more rapidly for the hotter 
gas). In contrast the fractional masses (Fig. 17b) in these temper- 
ature bands are negligible for |z| < 0.75 kpc, and increase above 
this height (less rapidly for the hotter gas). The temperature band 
5 X 10^ < T < 5 X 10^ K (green/black, dash-3dotted) is simi- 
larly distributed to the hotter gas (orange) in all profiles. It is how- 
ever identified with the warm phase, indicating that this is mainly 
hot gas cooling, a transitional state, which accounts for a relatively 
small volume fraction of the warm gas and especially a small mass 
fraction. The dramatic effect of increased resolution (Fig. 17a com- 
pared to c) is the significant increase in the very hot gas (red, dot- 
ted), particularly displacing the hotter gases (orange and green) but 
also to some degree the bulk warm gas (blue, dashed). This reflects 
the reduced cooling due to the better density contrasts resolved, as- 
sociating the hottest temperatures to the most diffuse gas. 

The middle temperature range 5x10^ <T<5xlO^K has 
a distinctive profile in both fractional volume and fractional mass, 
with minima near the mid-plane and maxima at about \z\ c:^ 400 pc, 
being replaced as the dominant component by hotter gas above 
this height. The fractional volume and vertical distribution of this 
gas is quite insensitive to the resolution. The distribution of the 
warm gas (5 x 10^ K ^ T < 5 x lO'^ K; blue, long-dashed) 
does not change much with increased resolution. However, the 
higher-resolution model has more of the cold phase (T < 500 K; 
black, solid and dash-dotted) and, especially, of the very hottest gas 
(T ^ 5 X 10^ K; red, dotted), at the expense of the intermediate 
temperature ranges. 

This can also be seen in the gas density and temperature prob- 



ability distributions shown in Fig. 18(a), (b): increased resolution 
modestly increases the abundance of cold gas and significantly en- 
hances the amount of very hot gas. The minima in the distributions 
(at density 10~^ cm~^, and at temperatures 10^ and 3 x 10^ K) 
appear independent of resolution, suggesting that the phase separa- 
tion is physical, rather than numerical. The distributions are most 
consistent in the thermally unstable range 313-6102 K. Higher res- 
olution also reduces the minimum further about the unstable range 
above 10^ K, as the highest temperature gas has lower losses to 
thermal conduction. The mean temperatures of the cold gas (60 K) 
warm gas (10^ K) and the mean warm gas density (0.14 cm~^) also 
appear to be independent of the resolution. However the natural log 
mean /Xn is about -8 for the hot gas, both within and without 2 pc of 
the mid-plane (with larger standard deviation for the gas near the 
mid-plane). This compares with values of -6.97 and -5.78 in our 
model with 4 pc resolution; i.e. a factor of about 1/3. This reflects 
the improved resolution of low density in the remnant interiors. 

The density and temperature probability distributions for 
WSWa are similar to those obtained by Joung & Mac Low (2006, 
their Fig. 7), who used a similar cooling function, despite the dif- 
ference in the numerical methods (adaptive mesh refinement down 
to 1.95 pc in their case). With slightly different implementation of 
the cooling and heating processes, again with adaptive mesh re- 
finement down to 1.25 pc, de Avillez & Breitschwerdt (2004, their 
Fig. 3) found significantly more cool, dense gas. It is noteworthy 
that the maximum densities and lowest temperatures obtained in 
our study with a non-adaptive grid are of the same order of mag- 
nitude as those from AMR-models where the local resolution is 
up to three times higher. At 4 pc our mean minimum temperature 
is 34 K, within the range 15-80 K for 0.625-2.5 pc (de Avillez & 
Breitschwerdt 2004, their Fig. 9). For mean maximum gas number, 
our 122 cm~^ is within their range 288-79 cm~^. 

The vertical density profiles obtained under the different nu- 
merical resolutions are shown in Fig. 19a. Although the density dis- 
tribution in Fig. 18a reveals higher density contrasts with increased 
resolution, there is little difference in the ^-profiles of the models. 
The mean gas number density at the mid-plane, n(0) — which with 
our course grid resolution excludes the contribution from Hll — is 
about 2.2 cm~^: double the observation estimates summarised in 
Ferriere (2001). This might be expected in the absence of the mag- 
netic and cosmic ray components of the ISM pressure, to help sup- 
port the gas against the gravitational force. 

However the vertical pressure distributions are consistent with 
the models of Boulares & Cox (1990, their Figs. 1 and 2), which 
include the weight of the ISM up to |z| = 5 kpc. The total pressure 
P(0) ^ 2.5 (2.0) X 10"^^ dyncm"^ for the standard (high) res- 
olution model is slightly above their estimate of about 1.9 for hot, 
turbulent gas. For the turbulent pressure alone we have po{0) ~ 
6.3 (7.9) X 10"^^ dyncm-2 falling to 1.0 (0.6) at \z\ = 500 pc 
and then remaining reasonably level. The pressures are generally 
slightly reduced with increased resolution, except for po near the 
mid-plane. Small scales are better resolved, so the turbulent struc- 
tures are a stronger component of the SN active region. These pres- 
sures are consistent with Boulares & Cox (1990), even though our 
model does not explicitly include the pressure contributions from 
the ISM above 1 kpc. 

Comparing our thermal pressure distribution (Fig. 18c) with 
de Avillez & Breitschwerdt (2004, their Fig. 4a) and Joung et al. 
(2009, their Fig. 2), the three models peak at 3.16, 1.3 and 4.1 x 
10~^^ dyncm"^, respectively. The latter models include \z\ — 
10 kpc and resolution up to 1.25 pc. Our data summarise the vol- 
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ume within z±l kpc, while the comparisons are within 10 kpc and 
125 pc respectively. 

We conclude that the main effects of the increased resolution 
are confined to the very hot interiors and to the thin shells of SN 
remnants; the interiors become hotter and the SN shell shocks be- 
come thinner with increased resolution (see Appendix B). Simulta- 
neously, the higher density of the shocked gas enhances cooling, 
producing more cold gas and reducing the total thermal energy. 
Otherwise, the overall structure of the diffuse gas is little affected: 
the probability distributions of thermal pressure are almost indistin- 
guishable, with our standard resolution fractionally higher pressure 
(Fig. 18c). 

We are satisfied that the numerical resolution of the reference 
model, A = 4 pc, is sufficient to model the diffuse gas phases re- 
liably. This choice of the working numerical resolution is further 
informed by tests involving the expansion of individual SN rem- 
nants (presented in Appendix B). 



9 DISCUSSION AND CONCLUSIONS 

The multi-phase gas structure obtained in our simulations appears 
to be robust, with overall parameters relatively insensitive to the 
physical (Section 4) and numerical (Section 8.3) details, includ- 
ing the parameterizations of the radiative cooling tested here (Sec- 
tion 8.1). We have identified natural temperature boundaries of the 
major phases using the variation, with height above the mid-plane, 
of the fractional volume occupied by the gas in relatively narrow 
temperature ranges. This confirms that the system can be satisfac- 
torily described in terms of just three major phases with temper- 
ature ranges T < 5 x 10^ K, 5 x 10^ ^ T < 5 x 10^ K and 
5 X 10^ ^ T < 5 X 10^ K. The most probable values of the 
variables we have explored (gas density, thermal and total pressure, 
perturbed velocity and Mach number) are practically independent 
of the cooling function chosen (Fig. 15). Moreover, this is true for 
the cold, warm and hot phases separately. A 3D rendering of a snap- 
shot of the density distribution from the reference model WSWa is 
illustrated in Fig 20, showing the typical location and density com- 
position of each phase separately. 

A conspicuous contribution to various diagnostics — espe- 
cially within 200 pc of the mid-plane, where most of the SNe are 
localised — comes from the very hot gas within SN remnants. Re- 
garding its contribution to integrated gas parameters, it should per- 
haps be considered as a separate phase. 

The fractional volume occupied by each phase is a convenient 
diagnostic and an important physical parameter. We have clarified 
the relation between the fractional volume and various probabilis- 
tic measures of a random distribution of density (or of any other 
quantity), and established an exact relation between the fractional 
volume and various density averages obtainable observationally (in 
Section 5). This represents a significant improvement upon the as- 
sumption of locally homogeneous gas, the only analytical tool used 
to date in determinations of the fractional volumes of the phases. 

The correlation scale of the random flows is obtained in Sec- 
tion 6, from the autocorrelation functions of the velocity compo- 
nents. Within 200 pc of the mid-plane, the horizontal velocity com- 
ponents have a consistent correlation scale of about 100 pc. In con- 
trast, the scale of the vertical velocity (which has a systematic part 
due to the galactic outflow of hot gas) grows from about 100 pc 
at the mid-plane to nearly 200 pc at z = 200 pc, and may do so 
further at larger heights (cf. Korpi et al. 1999). This is due to the 
increase of the fractional volume of the hot gas with distance from 



the mid-plane. Ai\z\ ~ 1 kpc most of the volume is occupied by 
the hot gas. As the interstellar gas flows out of the galactic disc into 
the halo, it must expand, and the scale of the expanding regions 
may be expected to become comparable to 1 kpc at ^ 1 kpc. It 
would be helpful to obtain estimates of the horizontal correlation of 
the flow above ±1 kpc, so that modelling of the galactic fountain 
might be adequately formulated. 

We find clear indication of cold gas falling back towards the 
mid-plane at speeds of a few km/s, hot gas involved in vigorous out- 
flow away from the mid-plane, and some warm gas entrained in this 
outflow (Section 7). The outflow speed of the hot gas increases up 
to 100 km within 100 pc of the mid-plane and then slowly de- 
creases. In contrast, the mean vertical velocity of the warm gas in- 
creases linearly with up to 20 km towards the upper bound- 
aries of our domain. 

Given that probability densities for gas temperature and num- 
ber density, calculated for individual phases, are clearly separated, 
the probability densities for both thermal and total pressure (the 
sum of thermal and turbulent) are not segregated at all. Despite its 
complex thermal and dynamical structure, the gas is in statistical 
pressure equilibrium. Since the SN-driven ISM is random in nature, 
both total and thermal pressure fluctuate strongly in both space and 
time (albeit with significantly smaller relative fluctuations than the 
gas density, temperature and perturbation velocity), so the pressure 
balance is also statistical in nature. These might appear to be ob- 
vious statements, since a statistical steady state (i.e., not involving 
systematic expansion or compression) must have such a pressure 
balance. Deviations from thermal pressure balance and observa- 
tions of significant regions of gas within the classically forbidden 
thermally unstable range (300 - 6000 K), which is also evident in 
our probability distributions, may lead to conclusions of an ISM 
comprising a broad thermodynamic continuum in pressure disequi- 
librium (Vazquez-Semadeni 2012, discussion on The controversy). 
The only systematic deviations from pressure balance are associ- 
ated with the systematic outflow of the hot gas (leading to lower 
pressures), and with the compression of the cold gas by shocks and 
other converging flows (leading to somewhat increased pressures). 
Even this can be further reconciled if we allow for the global verti- 
cal pressure gradient (cf. Fig. 7). It is evident that phases are locally 
in total pressure equilibrium. 

An important technical aspect of simulations of this kind is 
the minimum numerical resolution A required to capture the basic 
physics of the multi-phase ISM. We have shown that A = 4 pc is 
sufficient with the numerical methods employed here (Section 8.3). 
In addition to comparing results obtained for A = 4 pc and 2 pc 
with our own code, we have satisfied ourselves that our results 
are consistent with those obtained by other authors using adaptive 
mesh refinement with maximum resolutions of 2 pc and 1.25 pc. 

As with all other simulations of the SN-driven ISM, we em- 
ploy a host of numerical tools (such as shock-capturing diffusivity) 
to handle the extremely wide dynamical range (10^ < T < 10^ K 
and 10~^ < n < 10^ cm~^ in terms of gas temperature and num- 
ber density in our model) and widespread shocks characteristic of 
the multi-phase ISM driven by SNe. Their detailed description can 
be found in Section 2.4. We have carefully tested our numerical 
methods by reproducing, quite accurately, the Sedov-Taylor and 
snowplough analytical solutions for individual SN remnants (Ap- 
pendix B). 

The major elements of the ISM missing from the models pre- 
sented here are magnetic fields and cosmic rays. Analysis of the 
structure of the velocity field and its interaction with the magnetic 
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Figure 20. 3D snapshots, from model WSWa, of gas number density in (a) the cold gas, (b) the warm gas, and (c) the hot gas. In each plot regions that are 
clear (white space) contain gas belonging to another phase. The phases are separated at temperatures 500 K and 5 x 10^ K. The colour scale for logn is 
common to all three plots. 



field, effects of rotation, shear and SN rates will be the subject of a 
future paper. 
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APPENDIX A: NOTATION 

Table Al contains most of the symbols used in the text and their 
explanation, arranged alphabetically. 

APPENDIX B: EVOLUTION OF AN INDIVIDUAL 
SUPERNOVA REMNANT 

The thermal and kinetic energy supplied by SNe drives, directly or 
indirectly, all the processes discussed in this paper. It is therefore 
crucial that the model captures correctly the energy conversion in 
the SN remnants and its transformation into the thermal and kinetic 
energies of the interstellar gas. As discussed in Section 2.2, the size 
of the region where the SN energy is injected corresponds to the 



adiabatic (Sedov-Taylor) or the snowplough stage. Given the mul- 
titude of artificial numerical effects required to model the extreme 
conditions in the multi-phase ISM, it is important to verify that the 
basic physical effects are not affected, while sufficient numerical 
control of strong shocks, rapid radiative cooling, supersonic flows, 
etc., is properly ensured. Another important parameter to be chosen 
is the numerical resolution. 

Before starting the simulations of the multi-phase ISM re- 
ported in this paper, we have carefully confirmed that the model 
can reproduce, to sufficient accuracy, the known realistic analyt- 
ical solutions for the late stages of SN remnant expansion, until 
merger with the ISM. The minimum numerical resolution required 
to achieve this in our model is A = 4 pc. In this Appendix, we con- 
sider a single SN remnant, initialised as described in Section 2.2, 
that expands into a homogeneous environment. All the numerical 
elements of the model are in place, but here we use periodic bound- 
ary conditions in all dimensions. 

The parameters xi and vi are as applied in Model WSWa 
for A = 4 pc, but reduced here proportionally for A = 2 and 
Ipc. The constant C ^ 0.01 used in Eq. (10) to suppress cool- 
ing around shocks is unchanged. This may allow excess cooling at 
higher resolution, evident in the slightly reduced radii in Fig. Bl. 
For Model WSWah, xi and vi were just as in Model WSWa; for 
future reference, they should be appropriately adjusted, as should 
C to better optimise higher resolution performance. 

Bl The adiabatic and snowplough stages 

The Sedov-Taylor solution, 

is accurately reproduced with our code at the resolution A = 4 pc 
or higher. Here R is the remnant radius, ^sn the explosion energy, 
po the ambient gas density, and k ^ 2.026 for 7 = 5/3 (Ostriker 
& McKee 1988). 

Modelling even a single remnant becomes more challenging 
when radiative cooling becomes important. Here we compare nu- 
merical results with two analytic solutions for an SN remnant ex- 
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Figure Bl. The shell radius R of an SN remnant versus time, shown in (a) linear and (b) logarithmic scales; (c) the corresponding expansion speed R. Frame 
columns 1-3 are for different ambient gas densities, po x lO^'^ g cm~^ = 1.0, 0.1, 0.01 from left to right. Numerical results obtained under three numerical 
resolutions are shown: A = 4pc (black, solid), 2pc (green, dashed) and 1 pc (orange, dash-dotted). Dotted lines are for the standard snowplough solution 
(B2) (dark blue) and its modification by Cioffi et al. (1998) (light blue). The horizontal line in Panels (cl)-(c3) shows the sound speed in the ambient ISM. 



panding into a perfect, homogeneous, monatomic gas at rest. The 
standard momentum-conserving snowplough solution for a radia- 
tive SN remnant has the form 

ll/4 



R — Rq 



to) 



(B2) 



where Rq is the radius of the SN remnant at the time to of the 
transition from the adiabatic stage, and Rq is the shell expansion 
speed at to. The transition time is determined by Woltjer (1972) as 
that when half of the SN energy is lost to radiation; this happens 
when 

Ro = 230 kms- (^)''" ( tJ^ V ; (B3) 



IQSi erg 



the transitional expansion speed thus depends very weakly on pa- 
rameters. 

Cioffi et al. (1998) obtained numerical and analytical solutions 
for an expanding SN remnant with special attention to the transition 
from the Sedov-Taylor stage to the radiative stage. These authors 
adjusted an analytical solution for the pressure-driven snowplough 
stage to fit their numerical results to an accuracy of within 2% and 
5% in terms of R and R, respectively. (Their numerical resolution 
was 0.1 pc in the interstellar gas and 0.01 pc within ejecta.) They 
thus obtained 



R — Rr> 



4t_ 
3tp 



3/10 



(B4) 



where the subscript p denotes the radius and time for the transition 
to the pressure driven stage. The estimated time of this transition is 

3/14 



tp 



13 Myr 



( ) 



-4/7 



10^^ erg y vicm~ 
For ambient densities of po = (0.01,0.1,1) x 10~^^gcm~^, 



this yields transition times tp ^ (25, 6.6, 1.8) x 10^ yr and shell 
radii Rp ^ (130, 48, 18) pc, respectively, with speed Rp = 
(213,296,412) kms"^ 

This continues into the momentum driven stage with 



R_y _ 3.63 (t-tn 



1.29 



tp 



Rp 



(B5) 

where subscript m denotes the radius and time for this second tran- 
sition. 



V 10^1 erg 



-3/14 



/ no 

vlcm-37 



where i?ej — 5000 km is the initial velocity of the 4M0 
ejecta. For each po = (0.01,0.1,1.0) x 10~^^gcm"^, the 
transitions occur at tm = (168, 16.8, 1.68) Myr, and Rm = 
(1014, 281, 78) pc, respectively. The shell momentum in the lat- 
ter solution tends to a constant, and the solution thus converges 
with the momentum-conserving snowplough (B2); but, depending 
on the ambient density, the expansion may become subsonic and 
the remnant merge with the ISM before Eq. (B2) becomes applica- 
ble. 

We compare our results with the momentum-conserving 
snowplough solution and those of Cioffi et al. in Fig. Bl, testing 
our model with numerical resolutions A = 1, 2 and 4pc for the 
ambient gas densities po = (0.01,0.1,1.0,2.0) x 10"^^gcm~^. 
Shown in Fig. Bl are a linear plot of the remnant radius R ver- 
sus time to check if its magnitude is accurately reproduced, a dou- 
ble logarithmic plot of R{t) to confirm that the scaling is right, 
and variation of the expansion speed with time to help assess more 
delicate properties of the solution. We are satisfied to obtain good 
agreement with the analytical results for all the resolutions investi- 
gated when the ambient gas number density is below 1 cm~^. For 
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Table Al. Most important variables used in the text. 



Symbol Meaning 

Cs Adiabatic speed of sound 

Cp Heat capacity at constant pressure [ kpc^ Gyr~^ K~ ^] 

C Velocity autocorrelation function, Eq. (29) 

D/Dt Advective derivative, Eq. (5) 

V Velocity structure function, Eq. (28) 

eth Energy density, subscript thermal: 'th', kinetic: 'kin' 

Eq-^ Total energy injected into the ISM by a single SN 

fM,i Fractional mass of gas in phase i, Eq. (32) 

fY,i Fractional volume occupied by the phase i, Eq. (14) 

Qz Vertical acceleration due to the Galactic gravity, Eq. (7) 
Scale height of the Type I SN distribution. Section 2.2 

hii Scale height of the Type II SN distribution. Section 2.2 

k-Q Boltzmann's constant 

K Thermal conductivity (= Cppx) 

lo Velocity correlation scale 

mp Proton mass 

J\4 Mach number 

n Gas number density 

Hi Gas density averaged within a given phase i, Eq. (17) 

(rii) Gas density averaged over volume V of phase i, Eq. (18) 

p Thermal pressure 

P Total pressure (thermal plus turbulent) 

V Probability density 

rsN Characteristic radius of the SN energy injection site. Sec- 
tion 2.2 

rms Root-mean square 

s Specific entropy, [ergg~^K~^] 

Sn Parameter of the lognormal probability distribution, Eq. (13) 

S Velocity shear rate due to differential rotation 

SN Supernova (also as a subscript) 

T Gas temperature 

V Total volume of a region in Section 5.1 

Vi Volume occupied by an ISM phase labelled i. Section 5.1 

u Velocity perturbation: deviation of the gas velocity from the 

background rotational flow 

UQ Random velocity, p. 3 

U Large-scale shear flow (differential rotation), p. 3 

W Rate of strain tensor, Eq. (4) 

r Specific rate of photoelectric heating, [ erg g~ ^ s~ ^ ] 

A Numerical mesh separation (resolution of a simulation) 

Czy Shock-capturing viscosity 

Cx Shock-capturing thermal diffusivity 

A Radiative cooling rate, [ erg g~^ cm~^] 

)Li Molecular weight 

jjin Parameter of the lognormal probabihty distribution, Eq. (13) 

h> Kinematic viscosity 

Type I SN rate per unit surface area. Section 2.2 

^11 Type II SN rate per unit surface area. Section 2.2 

ft Angular velocity of the Galactic rotation 

(f)i Phase filling factor within the ISM phase i, Eq. (15) 

Volume filling factor of the ISM phase Eq. (16) 

PSN Rate of mass injection, per unit volume, by SNe, Section 2.2 

p Gas density 

&SN Rate of energy injection by SNe (per unit volume), as ki- 
netic energy in Eq. (2) and as thermal energy in Eq. (3), see 
Section 2.2 

af Variance of ISM phase i, Section 5 

Tcooi Radiative cooling time 

<^ Gravitational potential 

X Thermal diffusivity 



A = 4pc, the remnant radius is accurate to within about 3% for 
po = 10~^^ gcm~^ and underestimated by up to 6% for po = 
10~^^ g cm~^ . At higher numerical resolutions, the remnant radius 
is underestimated by up to 7% and 11% for po = 10~^^ gcm~^ 
and 10~^^gcm~^, respectively. For po = 10~^^gcm~^, excel- 
lent agreement is obtained for the higher resolutions, A = 1 and 
2 pc; simulations with A = 4 pc overestimate the remnant radius 
by about 20-25% in terms of and at t = 2 Myr. We empha- 
size that a typical SN explosion site in the models described in the 
main part of the paper has an ambient density no < 1 cm~^ so that 
A = 2, or 4 pc produce a satisfactory fit to the results, despite the 
much finer resolution of the simulations, of Cioffi et al. 

The higher than expected expansion speeds into dense gas can 
be explained by the artificial suppression of the radiative cooling 
within and near to the shock front as described by Eq. (10). Our 
model reproduces the low density explosions more accurately be- 
cause the shell density is lower, and radiative cooling is therefore 
less important. 

B2 The structure of the SN remnant 

Cuts through the simulated SN remnant are shown in Fig. B2 
for gas density, temperature and velocity, obtained for resolution 
A = 4pc and with ambient density po = 10~^^ gcm~^. In the 
temperature and velocity panels, we also include the profile of the 
shock viscosity from Eq. (9) (black dotted line), scaled to fit each 
plot. The temperature panels also show where net cooling is ap- 
plied to the remnant, T~^(r - pA) < from Eq.(3) (blue dashed 
line), while the velocity panels also show the ambient sound speed 
(pink dashed lines). The top panel depicts the initial distributions, 
at t = 0, with which the mass of 4Mq and 5 x 10^° erg each of 
thermal and kinetic energy are injected. The other panels are for 
t = 0.72 and 1.02 after the start of the evolution, from top to bot- 
tom, respectively; the actual simulation continued tot = 1.32 Myr, 
when the remnant radius reached 130 pc. 

The position of the peak of the density profile is used to de- 
termine the shell radius shown in Fig. Bl. The Rankine-Hugoniot 
jump conditions are not very well satisfied with the numerical pa- 
rameters used here. This is due to our numerical setup, essentially 
designed to control the shocks by spreading them sufficiently to be 
numerically resolvable in production runs that contain many inter- 
acting shocks and colliding SN shells. Better shock front profiles 
have been obtained with other choices of parameters and cooling 
control, and with better resolution. The density and temperature 
contrasts across the shock fronts are reduced by the shock smooth- 
ing, which inhibits the peak density and enhances gas density be- 
hind the shocks. In an isolated remnant, the peak gas number den- 
sity does not exceed 10cm~^, but in the full ISM simulation we 
obtain densities in excess of 100 cm~^, as a result of interacting 
remnants and highly supersonic flows. 

The interior of the SN remnant, if more dense due to numer- 
ical smoothing about the shock profile, would cool unrealistically 
rapidly, so that the SN energy would be lost to radiation rather than 
agitate the ambient ISM. The centre panels in Fig. B 1 clarify how 
the cooling suppression described in Eq. (10) reduces the cooling 
rate in the relatively homogeneous interior of the remnant, while 
still allowing rapid cooling in the dense shell where the gradient of 
the shock viscosity is small. It is evident from the temperature cuts 
that the remnant still contains substantial amounts of hot gas when 
its radius reaches 100 pc, so it would be merging with the ISM in 
the full simulation. 

The panels in the right column of Fig. Bl demonstrate that 
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Figure B2. One-dimensional cuts through the origin of an SN remnant expanding into gas of ambient density po = 10~^^ gcm~^, simulated with the 
numerical resolution A = 4 pc. The variables shown are (al)-(cl) gas number density (blue, solid), (a2)-(c2) temperature (red, solid), and (a3)-(c3) velocity 
(green, solid). The shock viscosity profile of Eq. (9) (scaled to fit the frame, black, dotted) is shown in the temperature and velocity panels; the net cooHng (blue, 
dashed), log(— (F — pA)+), from Eq. (3) is included in the temperature panel; and the ambient sound speed (pink, dotted) is also shown with the velocity. 
Panels in the top row (a) show the injection profiles used to initialise the remnant a.tt = 0; the lower panel rows are for the later times (h) t = 0.72 Myr and 
(c)t = 1.02 Myr. 



the interior gas velocity can be more than twice the shell speed. 
Due to the high interior temperature, this flow is subsonic, while 
the remnant shell expands supersonic ally with respect to its ambi- 
ent sound speed. The enhanced viscosity in the hotter interior (with 
viscosity proportional to the sound speed; see Section 2.4) inhibits 
numerical instabilities that could arise from the high velocities. In 
fact, accurate modelling of the SN interiors is not essential in the 
present context (where we are mainly interested in a realistic de- 
scription the multi-phase ISM), as long as the interaction of the 
remnant with the ambient gas is well described, in terms of the en- 
ergy conversion and transfer to the ISM, the scales and energy of 
turbulence, and the properties of the hot gas. 



APPENDIX C: BOUNDARY CONDITIONS AND 
NUMERICAL CONTROL OF ADVECTION AND 
DIFFUSION 

CI Top and bottom boundaries 

Unlike the horizontal boundaries of the computational domain, 
where periodic or sliding-periodic boundary conditions are ade- 
quate (within the constraints of the shearing box approximation), 



the boundary conditions at the top and bottom of the domain are 
more demanding. The vertical size of the galactic halo is of order of 
10 kpc, and nontrivial physical processes occur even at that height, 
especially when galactic wind and cosmic ray escape are impor- 
tant. As explained in Section 2.4, we do not attempt to model the 
full extent of the halo here. Therefore, it is important to formulate 
boundary conditions at the top and bottom of the domain that ad- 
mit the flow of matter and energy, while minimising any associated 
artifacts that might affect the interior. 

Stress-free, open vertical boundaries would seem to be the 
most appropriate, requiring that the horizontal stresses vanish, 
while gas density, entropy and vertical velocity have constant 
first derivatives on the top and bottom boundaries. These are im- 
plemented numerically using 'ghost' zones; i.e., three outer grid 
planes that allow derivatives at the boundary to be calculated in 
the same way as at interior grid points. The interior values of the 
variables are used to specify their ghost zone values. When a sharp 
structure approaches the boundary, the strong gradients are there- 
fore extrapolated into the ghost zones. This artificially enhances the 
prominence of such a structure, and may cause the code to crash. 
Here we describe how we have modified these boundary conditions 
to ensure the numerical stability of our model. 

To prevent artificial mass sources in the ghost zones, we im- 
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pose a weak negative gradient of gas density in the ghost zones. 
Thus, the density values are extrapolated to the ghost zones from 
the boundary point as 

y, ±Z ± kA) = (1 - A/0.1 kpc)y9(x, y, ±Z ± {k - 1) A) 

for all values of the horizontal coordinates x and y, where the 
boundary surfaces are Bi z = ±Z, and the ghost zones are at 
z = zbZ zb kA with /c = 1, 2, 3. The upper (lower) sign is used at 
the top (bottom) boundary. This ensures that gas density gradually 
declines in the ghost zones. 

To prevent a similar artificial enhancement of temperature 
spikes in the ghost zones, gas temperature there is kept equal to 
its value at the boundary, 

T(x, y, ±Z ± kA) = r(x, y, ±Z) , 

so that temperature is still free to fluctuate in response to the interior 
processes. This prescription is implemented in terms of entropy, 
given the density variation described above. 

Likewise, the vertical velocity in the ghost zones is kept equal 
to its boundary value if the latter is directed outwards, 

Uz{x,y,±Z ±kA) = Uz{x,y,±Z) , Uz{x,y,±Z) ^ 0. 

However, when gas cools rapidly near the boundary, pressure can 
decrease and gas would flow inwards away from the boundary. To 
avoid suppressing inward flows, where Uz{x, y, zbZ) ^ we use 
the following: if \uz{x, y, ±Z =F A)| < \uz{x, y, =bZ)|, we set 

Uz{x, y,±Z±A) = \ [uz{x, y, ±Z) + Uz{x, y, ±Z =F A)] ; 

otherwise, we set 

Uz{x, y, ±Z±A) = 2uz{x, y, ±Z) - Uz{x, y, ±Z ^ A) . 

In both cases, in the two outer ghost zones {k = 2, 3), we set 

Uz{x, y, ±Z ± kA) =2uz{x, y, ±Z ± (k - 1) A) 
-Uz{x,y,±Z±{k-2)A), 

so that the inward velocity in the ghost zones is always smaller 
than its boundary value. This permits gas flow across the boundary 
in both directions, but ensures that the flow is dominated by the 
interior dynamics, rather than by anything happening in the ghost 
zones. 

The Pencil code is non-conservative, so that gas mass is not 
necessarily conserved; this can be a problem due to extreme den- 
sity gradients developing with widespread strong shocks. Solving 
Eq. (1) for p, rather than In p, solves this problem for the snow- 
plough test cases described in Appendix Bl, with mass then being 
conserved within machine accuracy. However for the full model, 
once the ISM becomes highly turbulent, there remains some nu- 
merical mass loss. A comparison of mass loss through the vertical 
boundaries to the total mass loss in the volume indicates that nu- 
merical dissipation accounts for 1% per Gyr. The rate of phys- 
ical loss, from the net vertical outflow, was of order 15% per Gyr. 



C2 Time step control 

To achieve numerical stability with the explicit time stepping used, 
the CFL conditions have to be amply satisfied. For example, for 
advection terms, the numerical time step should be selected such 
that 



At < Hi 



A 



where Cs is the speed of sound, it = | it | is the amplitude of the 
perturbed velocity, i.e., the deviation from the imposed azimuthal 
shear flow C/, and At is a dimensionless number, determined empir- 
ically, which often must be significantly smaller than unity. Apart 
from the velocity field, other variables also affect the maximum 
time step, e.g., those associated with diffusion, cooling and heat- 
ing, so that the following inequalities also have to be satisfied: 



At < 



/tiA^ 



max(z/, 7x, ry) 



At < 



_K2_ 



where ki and K2 are further empirical constants and 



2^/|W|^ + C.(V-ti)^+Cx(V-^) 

cvT 



We use K, = Ki = 0.25 and K2 = 0.025. The latter, more stringent 
constraint has a surprisingly small impact on the typical time step, 
but a large positive effect on the numerical accuracy. Whilst the 
time step may occasionally decrease to below 0.1 or 0.01 years 
following an SN explosion, the typical time step is more than 100 
years. 

C3 Minimum diffusivity 

Numerical stability also requires that the Reynolds and Peclet num- 
bers defined at the resolution length A, as well as the Field length, 
are sufficiently small. These mesh Peclet and Reynolds numbers 
are defined as 



Pga 



uA 



ReA 



uA 



(CI) 



X X 

where Umax is the maximum perturbed velocity and A is the mesh 
length. For stability these must not exceed some value, typically 
between 1 and 10. Note that the Reynolds and Peclet numbers char- 
acterizing the flow are 25 times larger, since A = 0.004 is replaced 
by /o — 0.1 as the relevant turbulent length scale in the non-mesh 
quantities. 

In numerical modelling of systems with weak diffusivity, u 
and X are usually set constant, close to the smallest value consis- 
tent with the numerical stability requirements. This level strongly 
depends on the maximum velocity, and hence is related to the lo- 
cal sound speed, which can exceed 1500 km in our model. 
To avoid unnecessarily strong diffusion and heat conduction in the 
cold and warm phases, we scale the corresponding diffusivity with 
gas temperature, as T^^^. As a result, the diffusive smoothing is 
strongest in the hot phase (where it is most required). This may 
cause reduced velocity and temperature inhomogeneities within the 
hot gas, and may also reduce the temperature difference between 
the hot gas and the cooler phases. 

The effect of thermal instabihty is controlled by the Field 
length. 



A 



10~^^ erg cm^ 



where we have neglected any heating. To avoid unresolved den- 
sity and temperature structures produced by thermal instability, we 
require that Af > A, and so the minimum value of the thermal 
conductivity x follows as 



max(cs, u, U) ' 



1 



cool 



271 
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where TcooI is the minimum coohng time, and /3 is the relevant ex- 
ponent from the cooling function (e.g. as in Table 1 for WSW cool- 
ing). In the single remnant simulations of Appendix B, TcooI ^ 
0.75 Myr. In the full ISM simulations, minimum cooling times 
as low as 0.05 Myr were encountered. Xmin has maxima corre- 
sponding to /3 = 0.56, -0.2, -3, . . . for T = 313, 10^ 2.88 x 
10^ K, . . .. All of these, except for that at T = 313 K, result 
in Xmin < 4 X 10~^kms~^kpc at Cs = ci = lkms~^, so 
are satisfied by default for any xi sufficiently high to satisfy the 
PeA ^ 10 requirement. For T = 313 K, at Cs — c\ we have 
Xmin = 6.6 X 10~^kms~^ kpc > xi- Thus if cooling times as 
short as 0.05 Myr were to occur in the cold gas, we would have 
Af < A, and would be marginally under-resolved. Our analysis of 
the combined distribution of density and temperature, however, in- 
dicates that coohng times this short occur exclusively in the warm 
gas. 

With xi ~ 4.1 X 10~^ km kpc, as adopted in Section 2.4, 
then PeA ^ 10 is near the limit of numerical stability. (We dis- 
cuss our choice of thermal diffusivity further in Appendix D.) As 
a result, the code occasionally crashed (notably when hot gas was 
particularly abundant), and had to be restarted. When restarting, 
the position or timing of the next SN explosion was modified, so 
that the particularly troublesome SN that caused the problem was 
avoided. In extreme cases, it was necessary to increase x temporar- 
ily (for only a few hundred time steps), to reduce the value of PeA 
during the period most prone to instability, before the model could 
be continued with the normal parameter values. 



APPENDIX D: THERMAL INSTABILITY 

One of the two cooling functions employed in this paper, WSW, 
supports isobaric thermal instability in the temperature range 
313 < T < 6102 K where /3 < 1. (Otherwise, for the RBN cool- 
ing function or outside this temperature range for WSW cooling, 
we have /3 ^ 1 or F <C pA, so the gas is either thermally stable or 
has no unstable equilibrium.) 

Under realistic conditions of the ISM, thermal instability can 
produce very small, dense gas clouds which cannot be captured 
with the resolution A = 4 pc used here. Although the efficiency 
of thermal instability is questionable in the turbulent, magnetized 
ISM, where thermal pressure is just a part of the total pressure 
(Vazquez-Semadeni et al. 2000; Mac Low & Klessen 2004, and 
references therein), we prefer to suppress this instability in the 
model. However, we do that not by modifying the cooling func- 
tion, but rather by enhancing thermal diffusivity so as to avoid the 
growth of perturbations at wavelengths too short to be resolved by 
our grid. 

Following Field (1965), we introduce the characteristic wave 
numbers 



l)po>Cp 



7^CsTo 



- 1)Ct 

ncs 



kK 



IZcspo 



where 71 is the gas constant, and the derivatives Ct = (dC/dT)p 
and Cp = {dC/dp)T are calculated for constant p and T, respec- 
tively. The values of temperature and density in these equations. 
To and po, are those at thermal equilibrium, C{To, po) = with 
jC = pA — T. Isothermal and isochoric perturbations have the char- 
acteristic wave numbers kp and kr, respectively, whereas thermal 
conductivity K is characterised by kx- 

The control parameter of the instability is (p = kp/kx- 

The instability is suppressed by heat conduction, with the 



Table Dl. The unstable wavelengths of thermal instability, according to 
Field (1965), at thermally unstable equilibria (To , po) with the WSW cool- 
ing function. 
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largest unstable wave numbers given by (Field 1965) 
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for the condensation and wave modes, respectively, whereas the 
most unstable wave numbers are 



^mc — 
kuiw — 



p-1 



1/4 
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Table Dl contains the values of these quantities for the param- 
eters of the reference model WSWa, where we present the wave- 
lengths A = 2n/k rather than the wave numbers k. The unsta- 
ble wavelengths of thermal instability are comfortably resolved at 
To = 6102 K and 4000 K, with the maximum unstable wave- 
lengths Acc = 44 pc and 32 pc, respectively, being much larger 
than the grid spacing A = 4 pc. The shortest unstable wavelength 
of the condensation mode in our model, Acc = 5pcatT ^ 313K 
is marginally resolved at A = 4pc; gas at still lower tempera- 
tures is thermally stable. Unstable sound waves with Acw = 2 pc 
at T = 4000 K are shorter than the numerical resolution of the ref- 
erence model. However, for these wave modes to be unstable, the 
isentropic instability criterion must also be satisfied, which is not 
the case for /3 > 0, so these modes remain thermally stable. 

Thus, we are confident that the parameters of our models 
(most importantly, the thermal diffusivity) have been chosen so as 
to avoid any uncontrolled development of thermal instability, even 
when only the bulk thermal conductivity is accounted for. Since 
much of the cold gas, which is most unstable, has high Mach num- 
bers, thermal instability is further suppressed by the shock captur- 
ing diffusivity in the cold phase. 
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